幂次曲线

polyfit：最小二乘多项式曲线拟合

已知离散点上的数据集，即已知在点集上的函数值，构造一个解析函数（其图形为一曲线）使在原离散点上尽可能接近给定的值。

函数：p = polyfit(x,y,n)

其中：x是已知的离散数据点的横坐标，y是已知离散数据点的纵坐标，

n为需要拟合的最高次幂，由我们给定，运用不同的多项式进行拟合，

返回值p从左到右是高次到低次的多项式p(x)的系数，长度是n+1

p(x)=p1xn+p2xn−1+...+pnx+pn+1

函数：y=polyval(p,x); %根据拟合的函数得出x对应的因变量y的值

多项式n的阶数的确定：

可以用MATLAB的拟合工具箱，cftool进行选择

在MATLAB主窗口中输入 cftool 回车，会弹出拟合工具箱界面

选择拟合的参数，在右上角选择拟合方式为“Polynomial”，然后通过选择不同的degree，看右下角看离散点是否落在拟合曲线是以及Results里的SSE(方差)和R-square（相关系数），上图中的拟合，可以看到离散点都落在了拟合线上，并且相关系数为1，方差在-9的数量级上，拟合良好。

说明：大于一的幂级数与陡峭程度正相关，小于一的幂级数相当于曲线顺时针旋转90°，底数增减的量决定在横轴的平移。

<!DOCTYPE html>
<html lang="utf-8">
<meta http-equiv="Content-Type" content="text/html; charset=utf-8"/>
<head>
     <title>绘制二次三次曲线多次曲线</title>
    </head>

     <body οnlοad="draw()">
        <canvas id="myCanvus" width="1300px" height="640px" style="border:1px dashed black;">
            出现文字表示你的浏览器不支持HTML5
        </canvas>
     </body>
</html>
<script type="text/javascript">
<!--
    function draw(){
        var canvas=document.getElementById("myCanvus");
        var canvasWidth=1300;
        var canvasHeight=640;

        var context=canvas.getContext("2d");
        
        context.fillStyle = "white";
        context.fillRect(0, 0, canvasWidth, canvasHeight);

        context.strokeStyle = "black";
        context.fillStyle = "black";

        
        // 进行坐标变换：把原点放在左下角，东方为X轴正向，北方为Y轴正向
        var offsetY=320;// Y向偏移值，正值向上偏，用来画坐标轴
        var offsetX=650;// X向偏移值，正值向右偏，用来画坐标轴

        context.save();
        context.translate(0+offsetX,canvasHeight-offsetY);

        drawAxisXText(context);// 文字和线分开画比较好处理
        drawAxisYText(context);
        drawTitleText(context);

        context.rotate(getRad(180));
        context.scale(-1,1);        

        drawAxisX(context);        
        drawAxisY(context);       
        drawCurve(context);       

        context.restore();        
    }

    function drawTitleText(ctx){        
        ctx.lineWidth=0.5;
        ctx.strokeStyle='navy';
        ctx.fillStyle='navy';

        var x=350;
        var y=-250;

        // 写文字
        ctx.fillText("y=x^2 红色",x,y);    
        ctx.fillText("y=(x-3)^3 绿色",x,y+20);    
        ctx.fillText("y=(x-5)^4 黄色",x,y+40);    
        ctx.fillText("y=(x-7)^5 青柠色",x,y+60);    
        ctx.fillText("y=(x+3)^0.5 紫色",x,y+80);
        ctx.fillText("y=(x+5)^0.33 栗色",x,y+100);

        ctx.fillText("  绘制:逆火狂飙",x+170,y+30);
    }

    function drawCurve(ctx){
        var cds=[{}];
        var cds1=[{}];
        var cds2=[{}];
        var cds3=[{}];
        var cds4=[{}];
        var cds5=[{}];
        var cds6=[{}];

        var x,y,arr;
        for(x=-13;x<=13;x+=0.01){    
            y=Math.pow(x,2);//
            arr={"x":x,"y":y};
            cds.push(arr);
            
            y=Math.pow(x-3,3);//
            arr={"x":x,"y":y};
            cds1.push(arr);

            y=Math.pow(x-5,4);//
            arr={"x":x,"y":y};
            cds2.push(arr);

            y=Math.pow(x-7,5);//
            arr={"x":x,"y":y};
            cds3.push(arr);

            y=Math.pow(x+3,1/2);//
            arr={"x":x,"y":y};
            cds4.push(arr);

            y=Math.pow(x+5,1/3);//
            arr={"x":x,"y":y};
            cds5.push(arr);
        }

        paintCurve(ctx,"red",cds);
        paintCurve(ctx,"green",cds1);
        paintCurve(ctx,"yellow",cds2);
        paintCurve(ctx,"lime",cds3);
        paintCurve(ctx,"purple",cds4);
        paintCurve(ctx,"maroon",cds5);
        //paintCurve(ctx,"maroon",cds6);*/

        /*var ymax=-1000,ymin=1000,xmax,xmin;
        for(var i=0; i<cds.length; i++){  
            // 求y最大值
            if(cds[i].x<0 && cds[i].y>ymax){
                ymax=cds[i].y;
                xmax=cds[i].x;
            }

            // 求y最小值
            if(cds[i].x>=0 && cds[i].y<ymin){
                ymin=cds[i].y;
                xmin=cds[i].x;
            }
        } 

        console.log("ymin="+ymin+" xmin="+xmin+" ymax="+ymax+" ymin="+ymin+" xmax="+xmax);
        var SU=50;// Scale Unit
        // 极大值
        ctx.beginPath();
        ctx.moveTo(xmax*SU,ymax*5-5);
        ctx.lineTo(xmax*SU,ymax*5+5);

        ctx.save();
        ctx.scale(1,-1);
        ctx.fillText("ymax="+cutShort(ymax.toString(),8),xmax*SU,-ymax*5);
        ctx.restore();

        ctx.stroke();
        ctx.closePath();

        // 极小值
        ctx.beginPath();
        ctx.moveTo(xmin*SU,ymin*5-5);
        ctx.lineTo(xmin*SU,ymin*5+5);

        ctx.save();
        ctx.scale(1,-1);
        ctx.fillText("ymin="+ymin,xmin*SU,-ymin*5);
        ctx.restore();

        ctx.stroke();
        ctx.closePath();*/

        
    }

    function paintCurve(ctx,color,cds){
        var SU=50;// Scale Unit

        ctx.strokeStyle = color;
        ctx.beginPath();        
        for(var i=0; i<cds.length; i++){  
            ctx.lineTo(cds[i].x*SU,cds[i].y*SU);// 注意y轴比例
        }         
        ctx.stroke();
        ctx.closePath();
    }

    function drawAxisX(ctx){
        ctx.save();
        
        ctx.lineWidth=0.5;
        ctx.strokeStyle='navy';
        ctx.fillStyle='navy';

        var start=-650;
        var end=650;

        // 画轴
        ctx.beginPath();
        ctx.moveTo(start, 0);
        ctx.lineTo(end, 0);
        ctx.stroke();
        ctx.closePath();

        // 画箭头
        ctx.beginPath();
        ctx.moveTo(end-Math.cos(getRad(15))*10, Math.sin(getRad(15))*10);
        ctx.lineTo(end, 0);
        ctx.lineTo(end-Math.cos(getRad(15))*10, -Math.sin(getRad(15))*10);
        ctx.stroke();
        ctx.closePath();
        
        // 画刻度
        var x,y;
        y=5;
        for(x=start;x<end;x+=50){
            ctx.beginPath();
            ctx.moveTo(x, 0);
            ctx.lineTo(x, y);
            
            ctx.stroke();
            ctx.closePath();
        }

        ctx.restore();
    }

    function drawAxisXText(ctx){        
        ctx.lineWidth=0.5;
        ctx.strokeStyle='navy';
        ctx.fillStyle='navy';

        var start=-650;
        var end=650;

        // 写文字
        var x,y=5;
        for(x=start;x<end;x+=50){
            ctx.fillText(x/50,x,y+10);
        }
    }

    function drawAxisY(ctx){
        ctx.save();
        
        ctx.lineWidth=0.5;
        ctx.strokeStyle='navy';
        ctx.fillStyle='navy';

        var start=-300;
        var end=300;

        // 画轴
        ctx.beginPath();
        ctx.moveTo(0, start);
        ctx.lineTo(0, end);
        ctx.stroke();
        ctx.closePath();

        // 画箭头
        ctx.beginPath();
        ctx.moveTo(Math.sin(getRad(15))*10, end-Math.cos(getRad(15))*10);
        ctx.lineTo(0, end);
        ctx.lineTo(-Math.sin(getRad(15))*10, end-Math.cos(getRad(15))*10);
        ctx.stroke();
        ctx.closePath();
        
        // 画刻度
        var x,y;
        x=5;
        for(y=start;y<end;y+=50){// 注意y轴比例
            ctx.beginPath();
            ctx.moveTo(x, y);
            ctx.lineTo(0, y);
            
            ctx.stroke();
            ctx.closePath();
        }
    }

    function drawAxisYText(ctx){        
        ctx.lineWidth=0.5;
        ctx.strokeStyle='navy';
        ctx.fillStyle='navy';

        var start=-250;
        var end=350;

        // 写文字
        var x=-19,y=5;
        for(y=start;y<end;y+=50){

            if(y!=0){
                ctx.fillText(-y/50,x,y);// 注意y轴比例
            }
        }
    }

    function getRad(degree){
        return degree/180*Math.PI;
    }

    function cutShort(str,length){
        if(str.length>length){
            str=str.substr(0,length)+"...";
        }
        
        return str;
    }
//-->
</script>

本文转自张昺华-sky博客园博客，原文链接：http://www.cnblogs.com/xiandedanteng/p/8157169.html，如需转载请自行联系原作者

收集的问题：

如何用matlab来拟合幂律分布，怎样将拟合值和实际值进行对比，放在一个图中，又如何检验实际数据是否符合拟合函数。
如果不符合，如何来直接判断实际数据服从什么样的函数分布呢

在MATLAB里，对数据进行拟合，在双对数坐标下，看数据是否符合幂律分布，求出幂指数，并绘出图形。-In MATLAB, the pairs of data fitting, in double logarithmic coordinates to see whether the data meet the power-law distribution, find the power index, and draw the graph.以傅雷书信为例：（平移237个单位时间）

“>>”后是输入内容，行开头不含“>>”的是MATLAB运行的结果，“/”后是注释部分

>> x=[ ]; / 输入修正的书信间隔时间

>> y=[ ]; / 输入累积概率

>> loglog(x,y,'ko') / 画出双对数下坐标图（k表示黑色，o表示圆圈）

>> hold on / 保留刚才所画图表，以便继续在此图画出拟合直线

>> a=polyfit(log(x),log(y),1) / 求拟合直线的参数（一次项和常数项）

a =

-1.0700 5.9525

>> b=2.71828^5.9525 / 常数项转化

b =

384.7124

>> x=100:100000; / 根据上面所作的图指定x的取值范围

>> y=b*x.^ -1.0700 ; / 根据刚才的计算输入x和y的关系式

>> plot(x,y,'k') / 在双对数坐标下画出拟合直线

上面程序即可得到我们需要的图形，图形的再编辑可以在Figure窗口下的Edit-Figure Properties里修改（颜色、线条粗细、坐标轴命名等）。

但是得注意的是，用这个plotfit函数不太能够用来拟合很复杂的函数，而只是用来拟合线性的、二维之类的，而用它来拟合幂律分布的曲线时，只能考虑先截取一部分的点，然后用这个函数去拟合。

MATLAB软件提供了基本的曲线拟合函数的命令．

多项式函数拟合：P=polyfit(x,y,n)

其中n表示多项式的最高阶数，x，y为将要拟合的数据，它是用数组的方式输入．输出参数P为拟合多项式 P(1)*X^N + P(2)*X^(N-1) +...+ P(N)*X + P(N+1).的系数
多项式在x处的值y可用下面程序计算．
y=polyval(P,x,m)

线性：m=1, 二次：m=2, …

polyfit的输出是一个多项式系数的行向量。为了计算在xi数据点的多项式值，调用MATLAB的函数polyval。

例：

x=0:0.1:1;

y=[-0.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2];

A=polyfit(x,y,2)

Z=polyval(A,x);

Plot(x,y,’r*’,x,z,’b’)

polyfit不能保证你每次都能得到最优解，math的答案是使用数值计算。

个人认为，对于这种非线性的曲线，尽量不要使用ployfit, ployfit多项式抑合适合线性方程！！

用polyfit（）函数去拟合这么复杂的曲线不太合适，polyfit（）函数对于数据遵循多项式分布是比较好的，一般来说，利用polyfit（）函数拟合的阶数不要超过5阶。

如果是不需要得到拟合曲线的函数，只是把这些点利用一些光滑曲线连接，建议使用三次样条函数spline（）进行插值即可。

幂次变换的基本表达式为：y=cx ^r+b

其中c、r均为正数。与对数变换相同，幂次变换将部分灰度区域映射到更宽的区域中。当r=1时，幂次变换转变为线性变换。

(1) 当r<0时，变换函数曲线在正比函数上方。此时扩展低灰度级，压缩高灰度级，使图像变亮。这一点与对数变换十分相似。

(2) 当r>0时，变换函数曲线在正比函数下方。此时扩展高灰度级，压缩低灰度级，使图像变暗。

代码如下：

[cpp] view plain copy print ?

/******************************************************************************
* 作用: 幂次变换函数
* 参数: pDst 输出图像的像素数组
* 参数: pSrc 原始图像的像素数组
* 参数: nWidth 原始图像宽度
* 参数: nHeight 原始图像高度
* 参数: b 控制参数,表示曲线的上下偏移量
* 参数: c 控制参数,表示曲线的弯曲程度
* 参数: r 控制参数,表示函数的幂次
* 备注: 此函数对于彩色图同样适用
******************************************************************************/
int PowerTrans(BYTE* pDst, BYTE* pSrc, int nWidth, int nHeight, double b, double c, double r)
{
if (!pSrc || !pDst)
{
return EXIT_FAILURE;
}
// 映射表,用于256种灰度变换后的值
BYTE map[256];
// 保存运算后的临时值
double dTemp;
int i, j;
for (i = 0; i < 256; i++)
{
// 计算当前像素变换后的值
dTemp = c * pow(i / 255.0, r) * 255 + b;
// 如果超界则修改其值
if (dTemp < 0)
dTemp = 0.0;
else if (dTemp > 255)
dTemp = 255;
// 四舍五入
map[i] = int(dTemp + 0.5);
}
for (i = 0; i < nWidth * nHeight; i++)
{
for (j = 0; j < 4; j++)
pDst[i*4+j] = map[ pSrc[i*4+j] ];
}
return EXIT_SUCCESS;
}

WES幂曲线绘制实例

matlab幂函数曲线拟合代码-Geomech_IITR:用于标准井解释和地质力学建模的工具

幂函数型曲线裂纹平面问题的一般解 (2011年)

幂律分布拟合曲线

MATLAB polyfit曲线拟合及拟合最高次幂的选择

matlab幂函数曲线拟合代码-EchoDensity:参考实现和测试合成数据，用于声脉冲响应的排序中心时间回波密度...

绘制二次三次曲线多次曲线

论文研究-基于幂曲线映射的测量矩阵优化算法 .pdf

Matlab拟合曲线之幂律分布

点群6一维六方准晶中幂函数型曲线裂纹的反平面问题* (2012年)

若干含幂函数类对称曲线裂纹平面弹性问题的解析解 (2007年)

幂律分布的曲线方程拟合相关问题

matlab的双曲线代码-hp:双曲幂变换

非等间隔GM(1，1，tα)幂次时间项模型及其应用 (2015年)

图像处理-基本算法之幂次变换

Opencv 三次样条曲线(Cubic Spline)插值

l-曲线matlab代码-PowerLawsNanoscaleForces:纳米力幂定律的数据

python获取留存曲线的幂函数系数并计算生命周期LT

转：Matlab拟合曲线之幂律分布

数字图像处理(5)：幂次变换（C语言实现）

性感曲线幂指函数|x|^|x|超级应用

excel 方程拟合包含线性，指、幂函数，三阶贝塞尔绘图拟合曲线，贝塞尔（Bezier）求控制点

机器学习（幂次学院）笔记

任意次多项式曲线拟合（Delphi类）

c#chart已知一条曲线上的所有X和Y值，怎么求出这条曲线的幂

ISP（六）空间域图像变换（图像反转、对数变换、幂次变换、分段线性变换、直方图均衡与匹配）