matlab 小波包分解_[Cody精彩代码分享] 10的幂如何用4的幂来表达(MATLAB与Julia两个版本)...

Given an integer n, return the coefficients
c = [c_n,c_n-1,...,c_0]
Such that
10^n = c_n*4^(n) + c_n-1*4^(n-1) +...+ c_0*4^(0)
With the constraint that
c_n = c_0
For example,
n = 1
10^n = 10 = 2*4^(1) + 2*4^(0), so c = [2 2]

Expand 10^n to Powers of 4​ww2.mathworks.cn

问题分析:

从线性回归的角度来看, 这题的答案(系数)不唯一, 应该是有无数种答案.

因为回归自变量有n个(将4^n, 4^0合并为一个, 因为限制条件为c_n == c_0, 其他自变量为4^(n-1), 4^(n-2), ..., 4^1), 而样本个数只有一个(如果固定死n的话), 样本个数少于自变量个数, 属于典型的"欠定"问题, 理论上是有无数种答案的.

既然理论解上有无数种, 那么不妨构造出一种特殊的解.

既然样本个数为1, 那么理论上只需要一个自变量系数即可以拟合.

而限制条件正好要求c_n和c_0, 那么不妨让c_n和c_0(两者相等, 因此, 可以认为是一个自变量系数)作为要求的回归系数吧, 让其他系数等于0.

那么问题就简化成了:

10^n = c_n*4^n + c_0*4^0

= c_n*(4^n + 1)

=>

c_n = 10^n/(4^n + 1) (当 n >0时)

代码:

function ans = ten2pow4(n)
ans([1, n+1]) = 10^n/(4^n + sign(n));
end

注意里面的sign是为了应付特殊的情况, 即n == 0时,

当n == 0时,

10^0 = c_0*4^0

=>

c_0 = 10^0/4^0

比较前面的

c_n = 10^n/(4^n + 1)

发现区别无非在于分母多个"1"与少个"1"

使用sign, 就不需要用if语句来对两种情况做区分了.

等提交了代码, 发现更精彩的代码!!!:

size最小:

function ans = ten2pow4(n)
pow2(n)*diag(fliplr(pascal(n+1)));
end

size第二小:

function ans = ten2pow4(n)
arrayfun(@(i)nchoosek(n,i),0:n).*2.^n;

第一个代码第一眼不太好理解是什么意思. 经过我运行分析, 发现本质上是和第二个代码思路一样, 只是更简洁一些.

那么我就照着第二个代码讲解吧, 因为代码更直接, 更一目了然.

看看第二个代码, 我有些震惊, 为什么呢? 因为

1 代码里面没有任何"10", 也没有任何"4"(函数名除外, 哈哈), 但是竟然work了!

2 回归系数都是整数, 比起我的那个答案, 显然更具有美感.

经过了我的一些思考, 终于弄明白了原理, 并且写成文章分享给大家:

首先是那个"nchoosek(n,i)"和第一份代码中的"pascal", 这个很容易想到二项式定理.

二项式定理_百度百科​baike.baidu.com

(a + b)^n = sum(nchoosek(n, i) a^i*b^(n-i)), i为0到n.

那么就想办法将10分解成带有4的式子.

最容易想到的就是10 = (4 + 6), 很显然, leading solution里面的不是这样分解的, 当然, 这样分解也没有问题, 也是一种答案.

function ans = ten2pow4(n)
arrayfun(@(i)nchoosek(n,i).*(6.^(n - i)),n:-1:0);
end

但是, 无法满足约束条件: c_n == c_0.

理由是虽然nchoosek(n,i)是左右对称的, 但是6.^(n - i)显然不是左右对称的.

注意到代码里面带有"2^n",

容易想到了如何将10分解了:

10 = 2*(4 + 1)

=>

10^n = 2^n*(4 + 1)^n

对(4 + 1)^n进行二项式展开, 然后乘以常数2^n即可.

(4 + 1)与(4 + 6)的本质区别在于1的任何次方都是1, 显然是左右对称的, 从而能够保证c_n == c_0.

总结:

1 我的答案是从线性回归的角度出发, 构造出了大部分为0的系数.

2 leading solution的数学技巧更高, 答案更具有美感.

---2018-12-18--更新------------------------------------------------------------

最近学了点Julia, 现学现卖, 写一下Julia版的.

对

size第二小的MATLAB代码进行改写:

function ans = ten2pow4(n)
arrayfun(@(i)nchoosek(n,i),0:n).*2.^n;

改写后的Julia版:

ten2pow4(n) = binomial.(n, 0:n) .* 2^n

解释一下这行代码.

1 你没看错, 就一行, 连"function"都可以不用写. 因为一行代码的函数可以不用写"function", 这是Julia的一个语法糖, 让函数的定义更接近于数学语言.

2 binomial就是Julia版的nchoosek.

3 MATLAB通过使用arrayfun函数, 对不同内置向量化的函数(比如这里的nchoosek)进行向量化, 而Julia的方式不同, 通过函数名后面加".", 来实现向量化, 比如这里的"binomial.", 或者是运算符前面加".", 来实现向量化, 比如这里的".*", 这个方式和MATLAB类似.
4 注意Julia和MATLAB的一个重要区别: Julia的数字看起来是整数的话, 那么它就是真的是整数, 而MATLAB的很多时候, 看起来是整数, 实际上是浮点数(double型)

MATLAB默认都是double型的浮点数, 无论看起来是不是整数, 除非你显示定义它为整数型.

测试脚本我也改成了Julia版的:

ten2pow4(n) = binomial.(n, 0:n) .* 2^n

using LinearAlgebra
for n = 0:16
    coeff = ten2pow4(n);
    @assert(isequal(dot(coeff, 4 .^ reverse(0:n)),10^n) && isequal(coeff[1],coeff[end]) && coeff[1]>0)
    println(n, " passed.")
end

主要修改点:

1 将"assert"改为"@assert".

2 将索引的圆括号改为方括号.

吐槽一下MATLAB的圆括号太滥用了. 函数调用用圆括号, 索引也用圆括号, 有时候只看字面意思, 无法区分这两者.

3 "flip"改为了"reverse"

总结:

1 Julia除了速度快以外, 代码也可以写得很简洁. 当然这篇文章里面没有体现Julia的速度, 可以参考这篇文章:

菡父：[为什么学习Julia]Julia性能测试之“冰雹数”​zhuanlan.zhihu.com

2 Julia与MATLAB的语法有些接近, 你看, 我将MATLAB的运行脚本进行了少量的修改, 就可以让Julia运行了.

3 索引的方括号比MATLAb的圆括号更清晰, 消除歧义.

4 MATLAB的double型用得太多了, 经常还默认用来表示整数, 经常带来一些bug.

创作不易, 请大家"素质三连": 点赞, 收藏, 分享.