简介：FRED作为COM组件可以实现与Excel、VB、Matlab等调用来完成庞大的计算任务或画图，本文的目的是通过运行一个案例来实现与Matlab的相互调用，在此我们需要借助脚本来完成，此脚本为视为通用型脚本。&配置：在执行调用之前，我们需要在Matlab命令行窗口输入如下命令：enableservice('AutomationServer', true)enableservice('AutomationServer')结果输出为1，这种操作方式保证了当前的Matlab实体可以用于通信。&在winwrp界面，为增加和使用Matlab类型的目录库，我们需要如下步骤：1. 在FRED脚本编辑界面找到参考.2. 找到Matlab Automation Server Type Library3. 将名字改为MLAPP&&在Matlab里面有两种常用的数据发送选项PutWorkspaceData 及PutFullMatrix，PutWorkspaceData适用于存储一般的数据在工作区，并赋予其为变量，PutFullMatrix试用于复数数据。&图 编辑/参考&&现在将脚本代码公布如下，此脚本执行如下几个步骤：1. 创建Matlab服务器。2. 移动探测面对于前一聚焦面的位置。3. 在探测面追迹光线4. 在探测面计算照度5. 使用PutWorkspaceData发送照度数据到Matlab6. 使用PutFullMatrix发送标量场数据到Matlab中7. 用Matlab画出照度数据8. 在Matlab计算照度平均值9. 返回数据到FRED中&代码分享：&Option Explicit&Sub Main&& & Dim ana As T_ANALYSIS& & Dim move As T_OPERATION& & Dim Matlab As MLApp.MLApp& & Dim detNode As Long, detSurfNode As Long, anaSurfNode As Long& & Dim raysUsed As Long, nXpx As Long, nYpx As Long& & Dim irrad() As Double, imagData() As Double, reals() As Double, imags() As Double& & Dim z As Double, xMin As Double, xMax As Double, yMin As Double, yMax As Double& & Dim meanVal As Variant&& & Set Matlab = CreateObject(&Matlab.Application&)&& & ClearOutputWindow&& & 'Find the node numbers for the entities being used.& & detNode = FindFullName(&Geometry.Screen&)& & detSurfNode &= FindFullName(&Geometry.Screen.Surf 1&)& & anaSurfNode = FindFullName(&Analysis Surface(s).Analysis 1&)&& & 'Load the properties of the analysis surface being used.& & LoadAnalysis anaSurfNode, ana&& & 'Move the detector custom element to the desired z position.& & z = 50& & GetOperation detNode,1,move& & move.Type = &Shift&& & move.val3 = z& & SetOperation detNode,1,move& & Print &New screen position, z = & &z&& & 'Update the model and trace rays.& & EnableTextPrinting (False)& & & & Update& & & & DeleteRays& & & & TraceCreateDraw& & EnableTextPrinting (True)&& & 'Calculate the irradiance for rays on the detector surface.& & raysUsed &= Irradiance( detSurfNode, -1, ana, irrad )& & Print raysUsed & & rays were included in the irradiance calculation.&& & 'When using real number data to send to MATLAB, it is simplest to use PutWorkspaceData.& & Matlab.PutWorkspaceData(&irradiance_pwd&,&base&,irrad)&& & 'PutFullMatrix is more useful when actually having complex data such as with& & 'scalar wavefield, for example. Note that the scalarfield array in MATLAB& & 'is a complex valued array.& & raysUsed = ScalarField ( detSurfNode, -1, ana, reals, imags )& & Matlab.PutFullMatrix(&scalarfield&,&base&, reals, imags )& & Print raysUsed & & rays were included in the scalar field calculation.&&& & 'Calculate plot characteristics from the T_ANALYSIS structure. &This information is used& & 'to customize the plot figure.& & xMin = ana.posX+ana.AcellX*(ana.Amin-0.5)& & xMax = ana.posX+ana.AcellX*(ana.Amax+0.5)& & yMin = ana.posY+ana.BcellY*(ana.Bmin-0.5)& & yMax = ana.posY+ana.BcellY*(ana.Bmax+0.5)& & nXpx = ana.Amax-ana.Amin+1& & nYpx = ana.Bmax-ana.Bmin+1&& & 'Plot the data in Matlab with some parameters calculated from the T_ANALYSIS& & 'structure. &Set the axes labels, title, colorbar and plot view.& & Matlab.Execute( & surf(linspace(&&xMin &&,&&xMax &&,&&nXpx &&),linspace(&& yMin &&,& & yMax & &,& & nYpx & &),irradiance_pwd, 'EdgeColor', 'None');& )& & Matlab.Execute( &xlabel('X Position (& & GetUnits() & &)')& ) : Matlab.Execute( &ylabel('Y Position (& & GetUnits() & &)')& ) : Matlab.Execute( &zLabel( 'Irradiance' )& )& & Matlab.Execute( &title('Detector Irradiance')& )& & Matlab.Execute( &colorbar& )& & Matlab.Execute( &view(2)& )& & Print &&& & Print &Matlab figure plotted...&&& & 'Have Matlab calculate and return the mean value.& & Matlab.Execute( &irrad = mean(mean(irradiance_pwd));& )& & Matlab.GetWorkspaceData( &irrad&, &base&, meanVal )& & Print &The mean irradiance value calculated by Matlab is: & & meanVal&& & 'Release resources& & Set Matlab = Nothing&End Sub&最后在Matlab画图如下：&并在工作区保存了数据：&&并返回平均值：&与FRED中计算的照度图对比：& &例：&此例系统数据，可按照此数据建立模型&系统数据&&光源数据：Type: Laser Beam(Gaussian 00 mode)Beam size: 5；Grid size: 12；Sample pts: 100；相干光；波长0.5876微米，距离原点沿着Z轴负方向25mm。&对于执行代码，如果想保存图片，请在开始之前一定要执行如下代码：enableservice('AutomationServer', true)enableservice('AutomationServer')讯技光电2016年年度课程表：.cn/2016/VirtualLab Fusion建模方法剖析及光纤模拟 1月22日：.cn/2016/vf122.htmlFRED Optimum软件操作与实例剖析 3月21-22日 上海：.cn/2016/fred321.html扫一扫，关注讯技光电，了解更多软件信息！infotek(infotek)
文章为作者独立观点，不代表微头条立场
的最新文章
运用VirtualLab软件对于一个具有发散角的而激光激光器准直系统的性能研究。当提及模拟激光二极管时，FRED软件具有极大的灵活性。在这篇应用笔记中，将会描述简单到详细的激光光源模型。通过展示几个熟悉而创新的应用，如前房角镜、毛细管中的激光诱导荧光和人类皮肤模型，FRED和生物医疗产业的相关性能得到最好的表达。FRED结合人机界面(GUI)，可任意建构几何图形，并可直接由此接口中获得其对象外观，并拥有可满足此一精密设计需求的强大计算引擎之能力。在这个例子中，我们来看一个LED手电筒的简单示例。该案例主要说明了周期性结构内电磁场分布的严格计算。出于此目的，研究了两种不同的系统：具有狭缝的矩形光栅和具有三角形状的介质光栅(等腰三角形)。宽视场光学光谱仪是视觉受限的光学光谱仪，它是为第一代Thirty Meter Telescope 仪器而设计的。本文记录了成像模块配置中杂散光分析的进展。在项目的这一阶段杂散光分析的目标是提供预期的杂散光背景的基线评估。在FRED中共有四种方法设计柱透镜。在这个FRED模型中，侧入式LED智能手机显示上是一个虚拟原型。通过沿着波导加入渐变扩散片，可以实现均匀照明。基于物理光学的光学建模技术变得必不可少，其也是未来光学设计软件开发中顺理成章的一步。这就要求我们将光线追迹推广并将其与衍射建模技术联合起来。本示例主要介绍了1：6二元衍射分束器元件的严格参数优化所有的显示都需要一个机械装置来发送这个光穿过或不通过一个创建的文本或图形以便于显示信息。为控制光输出，会用到一个或几个光学仪器来控制所期望的显示信息的角和位置。对这种应用，FRED已经证实了相比传统成型方法它可节省30-50%的时间和成本。本文描述了如何在FRED中模拟空间滤波器，内容适用于所有的相干光束通过小孔的情况。本应用说明介绍了两种模拟LED的方法，强调了一些有用的分析工具。这个示例解释了如何导入测量光谱数据以在VirtualLab中构建一个光源。本文就一个具有Lambertian散射特性以及用户定义好的镜向反射系数的表面来进行演示。本文讨论了在安装FRED之后不同的喜好设定来优化GUI体验及光线追迹/分析特性。这个示例解释了如何导入测量光谱数据以在VirtualLab中构建一个光源。AZ眼睛模型达到了基于平均临床数据所决定的轴上和轴外像差等级。FRED文件包含了此眼睛模型及几个用来分析它的光源，并包括一个基于想要的屈光度来调整模型的内嵌脚本。FRED作为COM组件可以实现与Excel、VB、Matlab等调用来完成庞大的计算任务或画图，本文的目的是通过运行一个案例来实现与Matlab的相互调用，在此我们需要借助脚本来完成，此脚本为视为通用型脚本。本文介绍了一个脚本，即wollastonCreator.frs，根据输入到基本对话框中的用户技术参数来创建一个沃拉斯顿棱镜，且允许三个方向上具有不同的尺寸。FRED作为COM组件可以实现与Excel、VB、Matlab等调用来完成庞大的计算任务或画图，本文的目的是通过运行一个案例来实现与Matlab的相互调用，在此我们需要借助脚本来完成，此脚本为视为通用型脚本。说明 本文描述了如何使用一张位图（bitmap）图像作为光源反射体，这里我们要用到FRED详细光源结构。作为在光学工程中，高效的计算结果很难通过强力光线追迹来获得。使用辐射测量学技术，在很短的时间内，可以有效并准确地执行杂散光，照度均匀性和自发热辐射的计算来追迹必要数量的光线。懂得像差在光学系统中形成的原因，可以极大的帮助我们校正产生的这些像差，达到很好的成像质量。在新一代的VirtualLab 软件的名字里增加了”融合(Fusion)”旨在强调我们在该软件中结合了多种不同的建模技术。本文演示了如何模拟一个3反射镜5倍无焦望远镜，讨论了如何使用孔径选择基准抛物面的离轴部分来定义离轴抛物面。在建立该模型的过程中，使用一个脚本来追迹沿系统光轴的“主光线”，并且输出表面上光线入射的垂直位置说明 本文描述了如何使用一张位图（bitmap）图像作为光源反射体，这里我们要用到FRED详细光源结构。作为本文描述了从一个光学模型中导入胶合双透镜的一般过程，以Edmund光学45089消色差透镜作为一个范例（使用了Zemax文件）。平安夜到了，美妙的钟声又要敲响，让它敲响你2016年的幸运之门，带给你一整年的健康平安、幸福快乐、如意吉祥。VirtualLab提供了一个非常灵活的离轴抛物面反射镜元件，可用于各种应用中，例如望远镜仿真。如果一个光学系统要在不同参数或参数组合下进行多次模拟，参数运行文档则是一个强大的工具。VirtualLab提供各种结果的显示选项，包括动画。该示例中介绍了一个用于检测焦点位置的工具。这个示例演示了如何以一种期望的方式来配置满足要求的一维图表。我们希望通过VirtualLab Fusion可以帮助您享受高效的工作，我们将不断增加工具和技巧，来优化VirtualLab Fusion的用户体验。由于新的分布式计算功能已经添加到FRED中，所以，我们认为有必要将版本14.110作为此次开发周期中两个计划版本中的首个进行开发发布。用户可以使用属性浏览器获取工作文档的附加信息；这个示例介绍了属性浏览器的大体结构。这个示例对VirtualLab Fusion的用户界面的总体结构进行了一个基本的介绍。这个示例演示了如何使用快速访问工具栏用于常用的命令。这个示例描述了VirtualLab全局选项（Global options）对话框的基本要素。这个示例解释了全局选项对话框的性能（Performance）标签下的配置选项。介绍了如何利用VirtualLab软件中的光线追迹引进行参数优化。此案例模拟了光通过非傍轴球面透镜的传播和在焦点区域的矢量效应。由于新的分布式计算功能已经添加到FRED中，所以，我们认为有必要将版本14.110作为此次开发周期中两个计划版本中的首个进行开发发布。摘要：目前，FRED温度敏感性的评价可使用脚本语言实现。本文演示了一个双折射材料的折射率随温度变化而变化脚本摘要：目前，FRED温度敏感性的评价可使用脚本语言实现。本文演示了一个双折射材料的折射率随温度变化而变化脚本本文您将会学到如下内容：1.透镜基本参数输入；2.优化变量与评价函数设定；3.优化；4.照度分析；课程从光学薄膜的基本原理出发，逐渐深入，内容涵盖了光学薄膜从设计到制造所涉及的所有课题，并从项目管理和生产管理的角度对课程进行了完善。举办地点与时间：主办单位：讯技光电科技（上海）有限公司协办单位：北京理工大学、德国Lighttrans公司授本课程将以微纳光学理论基础及当前科研热门与VirtualLab在微纳光学的建模及设计能力相结合的形式进行。infotek讯技光电科技（上海）有限公司是一家专业的光电产业研发顾问服务公司，提供一流的专业光学仿真分析软件，并为客户提供全面的光电产业研发、制造整体解决方案。热门文章最新文章infotek讯技光电科技（上海）有限公司是一家专业的光电产业研发顾问服务公司，提供一流的专业光学仿真分析软件，并为客户提供全面的光电产业研发、制造整体解决方案。如何优雅地使用Matlab？
对Matlab软件本身以及语言一直提不起兴趣来，如果有机会我总是会用其他语言替代Matlab。但考虑到老师和同学，Matlab又是同他人交流的唯一语言，请教一下如何优雅地使用Matlab？我大概有行左右的Matlab经验，但由于写的代码都是一些重复性的简单工作，我对Matlab的认识还停留在翻本工具书调用个函数解决问题的初级阶段。对Matlab的语言的整体设计思路也一点都没有头绪。对其软件本身，对比Mathematica也觉得并不好用。请问下，如何深入的了解Matlab的设计/工作原理，并且优雅的使用它？
按投票排序
Matlab是个软件，作为用户的我们根本无从知道其工作原理，顶多就是看看官方手册介绍。所以我只说一下如何优雅的使用Matlab，仅仅是个人经验之谈，欢迎打脸1. 多借鉴其他语言特性，主要就是Fortran语言的模块化思想，还有C++中对象的概念；2. 想编写一个简洁优雅的程序要经过这么一种循环：看书－&学到的知识用到自己项目中－&读别人代码，找到闪光点－&修改自己代码－&再去读书……下面详细解释下1. 作为工科生难免会自己建立一个模型什么的，例如我们就经常鼓捣各种CFD模型。Matlab提供了一种几乎是积木似的语言，里面有各种现成的数学函数，但是想要把积木搭成城堡也需要一点学问不是？对我们来说这学问就是Fortran的模块化语言。下面再说就有点离题太远了，楼主有兴趣可以去找个现成的模型去看；2. 虽然最早的Matlab是用Fortran语言编写的，但是Matlab也可以有对象。（现在Fortran也有了，时髦啊！）在过去我从不用对象的时候，每次做项目的后处理都要针对各个项目重新写：提取结果文件，数据处理（插值排序等等），画图，输出图像或视频……每次都是那个累啊，终于有一天我看了本c++，里面提到了对象有代码通用性的好处，有了过去编程的基础再加上对象的思想，仅仅一个礼拜我就看完了Matlab官方的手册，后来我把每次后处理函数全部写成几个主要对象方法的形式（网格对象，计算结果对象，实测结果对象等等），虽然每次还是要改，但是比过去每次都重新编写要轻松太多了。而且我也发现“增加对象方法的数量”＋“让方法功能减少”就会使对象越来越灵活……当初真是应该多听几门计算机的编程语言课啊3. 为什么要多读书我就不说了，书籍是人类进步阶梯，要是谁说“我Matlab学的非常好但是从没看过Matlab的书”的话请 ，我保证不打死他……4. 其次还要多读别人代码，因为有的东西书里没告诉你，这时候要是你智商又不够悟不出来，所以参考一些大神还是很有必要的。举个栗子：初学者应该都认识reshape函数，我刚学的时候也觉得这函数老方便了，你看它能把矩阵化成向量，也能把向量化成矩阵，所以有时候一个插值函数里用不下十几次，而且还有的在循环里用（循环上千次啊）。后来在一个程序里看到，当需要一个矩阵M按列排列成向量时候，直接用M(:)的形式。虽然读程序可能会看走眼，不过不要忘了在上千次循环里这样可以有效提高程序运行速度，而且Matlab的矩阵在内存里就是按列排列的，可以说M(:)的读取速度基本上是最快的……所以要写好代码，读程序也很重要，这里推荐去
逛逛，里面又许程序提供下载。虽然许多仍不完美，但是肯定能找到提高自己编程地方。听说过这么一句话，程序＝数据结构 ＋ 算法，Matlab是一种能够很直接的把数学方法转化为程序的语言，所以能够很方便的验证你的算法是否正确，假如你只是想把它当作计算器来使，那么多用就可以了。上面只是大致讲下自己学习的经历，反正楼主也没问什么具体问题，所以不算偏题吧……废话这么多，求给来个赞
MATLAB不仅仅是一种语言, 而是一套面向科研工作的工具集合. 学习MATLAB是为了解决问题的, 不是为了学而学. 所以, 我认为首先要明确你的学习目的和应用环境, 要用MATLAB来解决哪些方面的问题. 掌握了MATLAB的基本语法和相关工具箱的使用后, 可以去MATLAB官网的
查找一些与你的工作相关的例子和源码, 会开阔你的思路. 网上一些开源的MATLAB工具箱也是很好的学习材料.
step1 .把你用过的函数全都edit一遍,看看怎么写的,edit为空的都是内核级的,不用管了极大改善你的代码风格step2 做个大算法解决个大问题,比如n=2000+的MTSP,你会发现你写的代码效率如此之低,就会绞尽脑汁的改善,自然就会查很多东西来学极大改善你的代码质量
edit看看那些内置函数的写法，可以用用vargarin之类的方法把自己写的函数改造成针对不同输入变量个数的函数，也可以学学GUI，为程序加上图形界面，感觉很有意思，也挺高端的
无所不能的Matlab！先上一张超有爱的图片~~M,你又调皮了^^=======================我是源码分割线 ========================[x,y]=meshgrid(-10:0.01:10);%% meshgrid是生成网格采样点的函数,该语句用于生成“心形”网格的采样点矩阵x,yz=-(17*x.^2-16*y.*abs(x)+17.*y.^2);%%该语句是
“心形” 函数，生成
“心形”网格的采样点矩阵 z[c,h]=contourf(z,100);%%contourf是绘制矩阵等高线的函数，该语句用于绘制 “心形”网格的等高线，100代表100条set(h,'linestyle','none')
完全向量化。可以用矩阵运算代替循环的时候，坚决不用循环。
优雅的使用？？学习python和r，回过头再来看matlab就会发现自己比他更优雅。。
我: 弄了一个算法完成模拟了我大学毕业前挂掉学分所有课程的仿真 感觉 优雅爆了
用MATLAB看时间
语法上Matlab和其他脚本语言没有太大不同，有覆盖各专业的科学计算库和丰富的数据结构。本身的GUI可以构建客户端程序，比直接开发PC客户端方便，并且Matlab也有UNIX版，可以在Mac和Linux上安装。楼主还是花一周时间啃本书学学吧，或许这是最好的结果。
已有帐号？
无法登录？
社交帐号登录苹果/安卓/wp
积分 350, 距离下一级还需 100 积分
权限: 自定义头衔, 签名中使用图片
道具: 彩虹炫, 雷达卡, 热点灯, 雷鸣之声, 涂鸦板, 金钱卡, 显身卡, 匿名卡下一级可获得
道具: 抢沙发
购买后可立即获得
权限: 隐身
道具: 金钱卡, 雷鸣之声, 彩虹炫, 雷达卡, 涂鸦板, 热点灯
各位童鞋：
& && &&&本人正在研究分位数回归在金融领域的应用，看了好多资料上的分位数回归都是采用R、SAS等其他统计软件进行分位数回归的参数估计，请问使用MATLAB如何进行分位数回归？MATLAB中有没有现成的分位数估计程序？对分位数回归感兴趣的同志可以加QQ一起讨论，谢谢！
载入中......
发表于10楼
function =quantreg(x,y,tau,order,Nboot)
% Quantile Regression
% USAGE: =quantreg(x,y,tau[,order,nboot]);
x,y: data that is fitted. (x and y should be columns)
Note: that if x is a matrix with several columns then multiple
linear regression is used and the "order" argument is not used.
tau: quantile used in regression.
order: polynomial orde ...
本帖最后由 zhangyiyiw 于
16:51 编辑
见过这个方法，可是自己不会用。有什么书籍介绍这个方法吗？
楼主自己看看把
10:32:48 上传
matlab分位数回归算法
投入天翼决图u币上的人都
liwenxue_137 发表于
楼主自己看看把写错了，这个只是分位数算法，稍微修改就可以了
投入天翼决图u币上的人都
liwenxue_137 发表于
楼主自己看看把谢谢你的M文件，不过这个M文件是求向量的分位数的，不是分位数回归，MATLAB应该没有现成的分位数回归程序。
没有查过，我用的是R软件做的，有个quantreg分位数包，应有尽有
投入天翼决图u币上的人都
我前一个月也用matlab，现在发现matlab做统计，很多包包都没有，就该R了。学习R只需要几天就可以了
投入天翼决图u币上的人都
liwenxue_137 发表于
我前一个月也用matlab，现在发现matlab做统计，很多包包都没有，就该R了。学习R只需要几天就可以了谢谢，stata也可以做分位数回归，我不是很熟悉R
function b = rq_fnm(X, y, p)
% Construct the dual problem of quantile regression
% Solve it with lp_fnm
[m n] = size(X);
u = ones(m, 1);
a = (1 - p) .*
b = -lp_fnm(X', -y', X' * a, u, a)';
function y = lp_fnm(A, c, b, u, x)
% Solve a linear program by the interior point method:
% min(c * u), s.t. A * x = b and 0 & x & u
% An initial feasible solution has to be provided as x
% Function lp_fnm of Daniel Morillo & Roger Koenker
% Translated from Ox to Matlab by Paul Eilers 1999
% Modified by Roger Koenker 2000--
% More changes by Paul Eilers 2004
% Set some constants
beta = 0.9995;
small = 1e-5;
max_it = 50;
[m n] = size(A);
% Generate inital feasible point
y = (A' \&&c')';
r = c - y * A;
r = r + 0.001 * (r == 0);& & % PE 2004
z = r .* (r & 0);
gap = c * x - y * b + w *
% Start iterations
while (gap) & small & it & max_it
& & it = it + 1;
& & %& &Compute affine step
& & q = 1 ./ (z' ./ x + w' ./ s);
& & r = z -
& & Q = spdiags(sqrt(q), 0, n, n);
& & AQ = A * Q;& && && & % PE 2004
& & rhs = Q * r';& && &&&% &
& & dy = (AQ' \ rhs)';& &% &
& & dx = q .* (dy * A - r)';
& & ds = -
& & dz = -z .* (1 + dx ./ x)';
& & dw = -w .* (1 + ds ./ s)';
& & % Compute maximum allowable step lengths
& & fx = bound(x, dx);
& & fs = bound(s, ds);
& & fw = bound(w, dw);
& & fz = bound(z, dz);
& & fp = min(fx, fs);
& & fd = min(fw, fz);
& & fp = min(min(beta * fp), 1);
& & fd = min(min(beta * fd), 1);
& & % If full step is feasible, take it. Otherwise modify it
& & if min(fp, fd) & 1
& && &&&% Update mu
& && &&&mu = z * x + w *
& && &&&g = (z + fd * dz) * (x + fp * dx) + (w + fd * dw) * (s + fp * ds);
& && &&&mu = mu * (g / mu) ^3 / ( 2* n);
& && &&&% Compute modified step
& && &&&dxdz = dx .* dz';
& && &&&dsdw = ds .* dw';
& && &&&xinv = 1 ./
& && &&&sinv = 1 ./
& && &&&xi = mu * (xinv - sinv);
& && &&&rhs = rhs + Q * (dxdz - dsdw - xi);
& && &&&dy = (AQ' \ rhs)';
& && &&&dx = q .* (A' * dy' + xi - r' -dxdz + dsdw);
& && &&&ds = -
& && &&&dz = mu * xinv' - z - xinv' .* z .* dx' - dxdz';
& && &&&dw = mu * sinv' - w - sinv' .* w .* ds' - dsdw';
& && &&&% Compute maximum allowable step lengths
& && &&&fx = bound(x, dx);
& && &&&fs = bound(s, ds);
& && &&&fw = bound(w, dw);
& && &&&fz = bound(z, dz);
& && &&&fp = min(fx, fs);
& && &&&fd = min(fw, fz);
& && &&&fp = min(min(beta * fp), 1);
& && &&&fd = min(min(beta * fd), 1);
& & % Take the step
& & x = x + fp *
& & s = s + fp *
& & y = y + fd *
& & w = w + fd *
& & z = z + fd *
& & gap = c * x - y * b + w *
& & %disp(gap);
function b = bound(x, dx)
% Fill vector with allowed step lengths
% Support function for lp_fnm
b = 1e20 + 0 *
f = find(dx & 0);
b(f) = -x(f) ./ dx(f);
这个就是matlab做分位数回归。最近在看这个，R，matlab都看了个够，2个月了
投入天翼决图u币上的人都
function [p,stats]=quantreg(x,y,tau,order,Nboot)
% Quantile Regression
% USAGE: [p,stats]=quantreg(x,y,tau[,order,nboot]);
%& &x,y: data that is fitted. (x and y should be columns)
%& && &&&Note: that if x is a matrix with several columns then multiple
%& && &&&linear regression is used and the &order& argument is not used.
%& &tau: quantile used in regression.
%& &order: polynomial order. (default=1)
%& && && & (negative orders are interpreted as zero intercept)
%& &nboot: number of bootstrap surrogates used in statistical inference.(default=200)
% stats is a structure with the following fields:
%& && &.pse:& & standard error on p. (not independent)
%& && &.pboot:&&the bootstrapped polynomial coefficients.
%& && &.yfitci: 95% confidence interval on &polyval(p,x)& or &x*p&
% [If no output arguments are specified then the code will attempt to make a default test figure
% for convenience, which may not be appropriate for your data (especially if x is not sorted).]
% Note: uses bootstrap on residuals for statistical inference. (see help bootstrp)
% check also:
% EXAMPLE:
% x=(1:1000)';
% y=randn(size(x)).*(1+x/300)+(x/300).^2;
% [p,stats]=quantreg(x,y,.9,2);
% plot(x,y,x,polyval(p,x),x,stats.yfitci,'k:')
% legend('data','2nd order 90th percentile fit','95% confidence interval','location','best')
% For references on the method check e.g. and refs therein:
%&&Copyright (C) 2008, Aslak Grinsted
%& &This software may be used, copied, or redistributed as long as it is not
%& &sold and this copyright notice is reproduced on each copy made.&&This
%& &routine is provided as is without any express or implied warranties
%& &whatsoever.
if nargin&3
& & error('Not enough input arguments.');
if nargin&4, order=[]; end
if nargin&5, Nboot=200; end
if (tau&=0)||(tau&=1),
& & error('the percentile (tau) must be between 0 and 1.')
if size(x,1)~=size(y,1)
& & error('length of x and y must be the same.');
if numel(y)~=size(y,1)
& & error('y must be a column vector.')
if size(x,2)==1
& & if isempty(order)
& && &&&order=1;
& & %Construct Vandermonde matrix.
& & if order&0
& && &&&x(:,order+1)=1;
& && &&&order=abs(order);
& & x(:,order)=x(:,1); %flipped LR instead of
& & for ii=order-1:-1:1
& && &&&x(:,ii)=x(:,order).*x(:,ii+1);
elseif isempty(order)
& & order=1; %used for default plotting
& & error('Can not use multi-column x and specify an order argument.');
pmean=x\y; %Start with an OLS regression
rho=@(r)sum(abs(r-(r&0).*r/tau));
p=fminsearch(@(p)rho(y-x*p),pmean);
if nargout==0
& & [xx,six]=sortrows(x(:,order));
& & plot(xx,y(six),'.',x(six,order),x(six,:)*p,'r.-')
& & legend('data',sprintf('quantreg-fit ptile=%.0f%%',tau*100),'location','best')
& & clear p
& & return
if nargout&1
& & %calculate confidence intervals using bootstrap on residuals
& & yfit=x*p;
& & resid=y-
& & stats.pboot=bootstrp(Nboot,@(bootr)fminsearch(@(p)rho(yfit+bootr-x*p),p)', resid);
& & stats.pse=std(stats.pboot);
& & qq=zeros(size(x,1),Nboot);
& & for ii=1:Nboot
& && &&&qq(:,ii)=x*stats.pboot(ii,:)';
& & stats.yfitci=prctile(qq',[2.5 97.5])';
% function ps=invtranspoly(p,kx)
% order=length(p);
% fact=cumprod(1:order);
% kx=cumprod(repmat(kx,1,order));
% for ii=2:order
%& &&&for jj=1:ii-1
%& && && &N=order-jj+1;
%& && && &k=ii-
%& && && &k=min(k,N-k);
%& && && &ps(ii)=ps(ii)+p(jj)*kk(k)*fact(N)/(fact(k)*fact(N-k));
论坛好贴推荐
  |
  |
  |
  |
  |
  |
如有投资本站或合作意向，请联系（010-）；
邮箱：service@pinggu.org
投诉或不良信息处理：（010-）
京ICP证090565号
京公网安备号
论坛法律顾问：王进律师
