积分曲线的样条估计及在脑图像技术中的应用.doc

苏州大学本科生毕业设计 论文 目目 录录 摘要摘要 1 ABSTRACTABSTRACT 2 第第 1 1 章章 前言前言 3 1 1 脑图像技术的背景知识介绍 3 1 2 本文所涉及理论知识的简单介绍 3 1 3 本文的主要工作 4 第第 2 2 章章 理论知识详细介绍理论知识详细介绍 4 2 1 B 样条函数回归理论 5 2 2 欧拉方法 6 第第 3 3 章章 问题简述问题简述 8 第第 4 4 章章 模型的建立以及程序的编写模型的建立以及程序的编写 8 4 1 已知的积分曲线的选择 8 4 2 积分曲线最优估计的求解 10 4 2 1 B 样条基函数矩阵 M 文件的建立 10 4 2 2 的求解 10 x m 4 2 3 的求解 13 tx 4 3 协方差矩阵及椭圆置信区域的求解 14 4 3 1 协方差矩阵的求解 14 tG 4 3 2 椭圆置信区域的求解 15 第第 5 5 章章 结果分析结果分析 17 5 1 对积分曲线整体估计好坏的判断分析 17 5 2 基于置信度概念的估计好坏的判断 19 参参 考考 文文 献献 23 致致 谢谢 23 附附 录录 24 苏州大学本科生毕业设计 论文 1 积分曲线的样条估计及在脑图像技术中的应用积分曲线的样条估计及在脑图像技术中的应用 对已知积分曲线的模拟实验对已知积分曲线的模拟实验 摘要摘要 现代神经学中 神经科学家在研究神经纤维的位置时 常用的脑图像技术称为扩散张量 成像 神经纤维的位置是通过积分曲线建立模型 但是该积分曲线并不是直接观察到的 比 如说一个二维向量场是在一个被扰动的规则网格上观测到的 而规则网格是通过添加随机噪 声而被扰动的 关键的步骤是求出满足某常微分方程的积分曲线 但该微分方程的函数形式 通常是未知的 本文首先运用样条函数回归理论 通过光滑化从随机场数据中提取出微分方 程的函数形式 进而采用欧拉方法求出积分曲线的估计 再构造出相应的椭圆形置信区域 为了判断这种方法的好坏 本文选取了已知的积分曲线 对微分方程 2 1 sin 8 3 2 1 cos 8 3 2 1 ttx ttx 2 1 8 7 0 0 xxxm x t dt td 21 112 12 12 212 0 5 0 5 T T ij ij x tem x tx t x tx t x tem x tx t m 其中是的独立同分布的有界随机向量 进行相同的方法编程计 12 1 n ijijij i j ee e 0 ij E e 算 我们通过欧拉方法在点开始迭代 计算出该已知积分曲线的估计 解 2 1 8 7 0 0 xx 出在光滑向量场中的微分方程 同时 我们求出积分曲线的估计的协方差矩阵 来求得沿着 积分曲线的置信椭圆带 并在三维图中画出 最后 基于置信度的概念 对该方法程序进行 好坏分析 使得置信椭圆带较好地包括住已知的积分曲线 并且程序运行速度较快 关键词 关键词 B 样条 欧拉方法 置信度 置信椭圆区域 积分曲线 苏州大学本科生毕业设计 论文 2 ABSTRACTABSTRACT Mordern neuroscientists study the location of neural fibers the main brain imaging technique is called diffusion tensor imaging They are modeled by integral curves that are not observed directly Instead for example a two dimensional vector field is observed on a regular grid perturbed by additive random noise Therefore the pivotal step is to compute an integral curve that satisfies an ordinary differential equation ODE But the function form of the ODE is frequently unknown Firstly the article apply the spline function regression theory Through smoothing message from random data one can get the ODE s function form in order to calculate the estimator of the integral curve by Euler s method then construct the corresponding confidence ellipses In order to analyse advantages and disadvantages of the method choose a simulative integral curve We use the same method to calculate the 2 1 sin 8 3 2 1 cos 8 3 2 1 ttx ttx ODE 2 1 8 7 0 0 xxxm x t dt td 21 112 12 12 212 0 5 0 5 T T ij ij x tem x tx t x tx t x tem x tx t m where are i i d bounded random vectors with by 12 1 n ijijij i j ee e 0 ij E e programming We compute our estimator of the true integral curve starting at a given point by using Euler s method to numerically solve the ODE 2 1 8 7 0 0 xx generated by the smoothed field Simultaneously we track the covariance matrix of the estimator of the integral curve in order to obtain confidence ellipses along the integral curve which are then plotted in three dimensional space Lastly we analyse the method based on the concept of degree of confidence and improve the program in order to make the confidence ellipses tube contain the true integral curve well and the program run faster Key Words B spline Euler s method the degree of confidence confidence ellipses integral curves 苏州大学本科生毕业设计 论文 3 第第 1 1 章章 前言前言 1 11 1 脑图像技术背景知识介绍脑图像技术背景知识介绍 脑功能成像是使用脑局部成像技术 完全无损伤地通过图像扫描显示出脑的高级活动的 一种新方法 轴突 Axon 是由神经元 即神经细胞的细胞本体长出突起的 其功能是传 递细胞本体的动作电位至突触 神经系统中 轴突是主要的神经信号传递渠道 大量轴突牵 连一起 以其外型类似而称为神经纤维 现代神经学家在研究神经纤维的位置时 常用的脑 图像技术是扩散张量成像 在扩散张量成像中 表达器官对于水的在各个方向的微分透性的 张量可以用来产生大脑的扫描图 在脑图像技术中关键的一步是求出满足某常微分方程的积 分曲线 常用方法是根据非线性方程的牛顿迭代法构造一个用于非线性常微分方程初值问题 数值解的迭代解法 在2007年 Institute of Mathematical Statistics的研究人员Koltchinskii Sakhanenko and Cai主要运用局部光滑的方法从随机场数据中提取出微分方程的函数形式 主要是在每 个数据所对应点的小领域内来进行线性回归 虽然能较好的拟合出积分曲线的精确解 但是 随着随机数据量的不断增大 需要进行线性回归的次数也越来越多 从而导致计算机工作量 的大幅度增加 运行时间较长 不能高效地解决实际问题 所以本文将通过研究 设计出一个新的解法 运用计算机软件编写出程序 且不断的改 善 最后得到运算速度快且便于非数学类专家使用的方法 利于促进脑科学的研究 相信这 将给脑神经医学研究提供有力的统计工具 并且在核磁共振成像等物理方面也有重要的统计 运用 1 21 2 本文所涉及理论知识的简单介绍本文所涉及理论知识的简单介绍 插值方法插值方法 在科学研究和生产时间活动中所遇到的大量函数 有相当一部分是通过观测或实验得到 的 虽然其函数关系在某个区间 a b 上是存在的 但却不知道具体的解析表达式 xfy 只能通过观测或实验得到一些离散点的函数值 导数值 因此希望对这样的函数用一个简单 的解析表达式近似的给出整体上的描述 在数学上 即在某个区间 a b 上存在且连 xfy 苏州大学本科生毕业设计 论文 4 续 但不知道解析表达式 只给出 a b 上离散点的函数值 需要计算在 i x i xfy xfy a b 上其他点的函数值 完成这项工作的方法有几种 如拉格朗日插值 分段 3 次Hermite插值和三次样条插值 等 实践表明 拉格朗日插值和Hermite插值函数对于结点较多且拟合线型复杂的情况 只 能保证曲线的连续性 却不能保证曲线的光滑性 即各段连续点处导数的连续性不能保证 要保证曲线连续且曲率也连续 这就要求插值函数具有连续的二阶导数 对于基点较多且分 布不规则的神经纤维情况 做一个高次插值多项式是不理想的 因为它带有近似性 且计算 也相当复杂 这样 很自然地想到样条插值 由于 3 次B样条函数是多项式样条函数中的一 类具有正性与局部支撑性的样条函数 利用B样条函数来表示插值样条将大大地化简计算 而且计算的数值稳定性也很好 常微分方程初值解问题常微分方程初值解问题 微分方程是用数学方法分析客观现象变化规律的有力工具 建立微分方程只是解决实际 问题的第一步 通常需要求出方程的解来分析实际现象 并加以检验 虽然有时可以利用微 积分方法求出某些微分方程的解析解 但是实际上大量的微分方程都不能获得它们的解析表 达式 即使有时能获得解析解 其表达式也非常复杂 难以讨论其性质 因此必须通过数值 求解的方法算出微分方程在某些离散点出的近似解 进而分析微分方程所反映的客观规律 最常用的两种数值算法是欧拉 Euler 方法和龙格 库塔 Runge Kutta 方法 1 31 3 本文的主要工作本文的主要工作 本文主要是运用样条函数回归理论 通过光滑化从随机场数据中提取出微分方程的函数 形式 进而构造积分曲线的估计和相应的椭圆形置信区域 通过本文所研究出的解法 我们 需要快速 高效地得到微分方程函数形式 本文对已知的积分曲线 采用相同的方法 利用 MATLAB 编写出相应的程序 通过此 程序得到微分方程的函数形式 且能够在三维的空间里构造出已知积分曲线的最优估计和相 应的椭圆形置信区域 最后得到的三维图中不仅要求积分曲线的估计值和真实值的差趋于零 还要求置信椭圆带按指定的概率覆盖住真实曲线 第第 2 章章 理论知识的详细介绍理论知识的详细介绍 苏州大学本科生毕业设计 论文 5 2 12 1 B B 样条函数回归理论样条函数回归理论 定义定义 3 1 13 1 1 m m 次样条函数次样条函数 设 a b 上给定一个划分 bxxxa n 10 如果函数满足下列条件 xS 1 函数 1 baCxS m 2 函数在每个子区间上是 m 次代数多项式 xS 1 ii xx 2 1 0 ni 则称是关于结点划分的m次样条函数 xS 设满足上述条件的m次样条函数的全体组成的集合为 容易验证 是一 mS mS 个m n维的线性空间 定义定义 3 1 23 1 2 截断幂函数截断幂函数 定义截断幂函数为 0 0 0 x xx x m m 定理定理 3 1 13 1 1 中的n m个样条函数 mS 1 2 1 2 1 0 nixx mkx m i k 在区间 a b 上线性无关 从而可得出的维数为n m 即是n m维线性空间 mS mS 为了构造B 样条函数 先对划分加入新点扩展为 mnnnm xxbxxxaxx 1101 称向量为结点向量 101mnm xxxxx 令 m m xtxt 视为参数 是 的函数 当x xt m t mnnm xxxxxxt 101 时 都是关于划分的样条函数 记为 01 xxxxxxxx mnmmmmm xtt mm 关于所作的 m 1 阶差商记为 11 mjjj xxxt 11 mjjjm xxx 定义定义 3 1 33 1 3 m m 次次 B B 样条样条 苏州大学本科生毕业设计 论文 6 设是结点序列 令函数关于 i x m m xtt 1 txx mjmj 的 m 1 阶差商 11 mjjj xxxt 111 mjjjmjmjmj xxxxxB 1 1 nmmj 称为第个 m 次 B 样条函数 简称 B 样条函数 j 利用差商的性质 1 11 mj jk kjm m k mjjjm x xx xxx 其中 得 11 mjjjjm xtxtxtt 1 1 1 mj jk kjm m k jmjmj x xx xxxB 由 1 式定义的 n m 个样条函数是线性无关的 所以组成的一组基 这 xB mj mS 样对于定义于区间 a b 关于划分的 m 次样条函数都可以表示为 mSxS 2 1 n mj mjj xBaxS 这样 求划分上的样条函数的问题 就归结为求表达式 2 中的系数 xS 的问题 求系数一般归结为解线性方程组 11 nmm aaa j a 由 B 样条函数的定义 注意到差商的性质 可以得到 B 样条具有下列重要性质 1 递推关系 其他 0 1 1 0 jj j xxx xB mkxB xx xx xB xx xx xB kj jkj kj kj jkj j nj 2 1 1 1 11 1 1 2 B 样条函数的导数 当m 0时 当 m 1时除为结点外 时 0 0 xBj1 mx 11 1 11 jmj mj jmj mj mj xx xB xx xB mxB 2 22 2 欧拉方法欧拉方法 设有 3 00 yxy yxfxy 苏州大学本科生毕业设计 论文 7 其中已知函数对满足Lipschitz条件 即存在常数使 yxfyL 2121 yyLyxfyxf 以保证微分方程 3 的解存在且惟一 在满足 Lipschitz 条件下求解微分方程 3 xyy 称为一阶常微分方程的初值问题 所谓求微分方程 3 的数值解 就是计算 精确 解在一系列离散点 xy 上的近似值 通常选取相等的计算步长h 于是 n xxxx 210 1 0 nhxxn 2 1 n 欧拉方法的基本思想是在小区间上用数值微分的前差公式 1 nn xx 代替方程 1 中左端的导数 而方程右端函数中的取 h afhaf af y xyxfx 中的某一点 于是 1 nn xx 4 11 nnnn xxxxyxhfxyxy 取内不同的点 就得到以下不同的公式 x 1 nn xx 若中的取区间的左端点 即 将的近似值记作 xyxfx 1 nn xx n x nn xyxf n xy 即 则得到近似计算公式 n y nn xyy 11 nn xyy 5 2 1 0 1 nyxhfyy nnnn 称为向前欧拉公式 意思是从 算出 从 算出 依次向前一步步地计算 0 x 0 y 1 y 1 x 1 y 2 y 若中的取区间的右端点 类似地可得 xyxfx 1 nn xx 1 n x 6 2 1 0 111 nyxhfyy nnnn 称为向后欧拉公式 形式上是从 向后算出 但是初值问题需要从和计算 1 n x 1 n y n y n y 1 n x 当函数对非线性时 6 无法直接进行计算 通常要用迭代法求解 比如先 1 n y yxfy 有向前欧拉公式算出的初值 1 n y 0 1nnnn yxhfyy 再按如下迭代公式计算 2 1 0 11 1 1 kyxhfyy k nnn k n 若迭代过程收敛 则所要求的值 每一步都要进行迭代 11 lim k n k n yy 2 1 0 n 将 5 式 6 式平均 得到 7 2 1 0 2 111 nyxfyxf h yy nnnnnn 称为梯形公式 这里取了区间左右端点的函数的平均值 xyxf 1 nn xx 显然 7 式也是隐式公式 迭代过程计算量大 通常采用预测 校正方法加以改进 先用 5 式预测 再把预测值代入 7 式右端做一次校正 即 1 n y 苏州大学本科生毕业设计 论文 8 2 1 0 2 111 1 nyxfyxf h yy yxhfyy nnnnnn nnnn 称为改进欧拉公式 第第 3 3 章章 问题简述问题简述 现代神经学研究神经纤维的位置 常用的脑图像技术称为扩散张量成像 首先我们已知 一个映射 其中是一个有界开集 通过扩散张量 核磁共振成像 DT 2 Rm 2 R MRI 我们得到在一些带有随机误差的离散点处的映射值 不失一般性 我们取 2 1 0 并将观察到的数据表示成 12 1 ijijijijij YYi jn Ym xe 其中是的独立同分布的有界随机向量 为中的 12 1 n ijijij i j ee e 0 ij E e 12 ijijij xx x 2 1 0 网格点即 12 1 1 ijijij xxinjn x 科学终端用户主要感兴趣的是追踪二维的曲线即轴突的轨迹 0 xT 映射满足常微分方程 0 xT 0 0 0 dt ttT dt x m xxx 其中 均为已知的固定值 微分方程写成积分形式即为 0 x0 T 120 0 0 t tx tx tt ds tT xxm x 我们需要根据上面的方程以及估计出的来解出曲线 在已知数据m tx 估计出是量化大脑中的神经纤维带的轨迹的主要方法 nji ijij 1 Yx 0 t tT x DT MRI 数据的扰动会使神经纤维带的轨迹不清晰 而从业者需要理解这种扰动干扰轨迹的 程度 所以 除了对的估计 对的不确定性研究也很重要 tx tx 我们用的插件式估计使得 tx tt n XX 120 0 0 t tX tXts dstT Xxm X 其中 是的一个估计 估计好坏的理论评估是根据渐近的协方差来确定的 而此协方差 m m 是可以通过估计曲线的逐点置信椭圆来判断 tX 第第 4 4 章章 模型的建立以及程序的编写模型的建立以及程序的编写 4 14 1 已知的积分曲线的选择已知的积分曲线的选择 苏州大学本科生毕业设计 论文 9 本文主要采用已知的积分曲线 利用 MATLAB 编写出相应的程序 通过此程序得到微 分方程的函数形式 且能够在三维的空间里构造出积分曲线的最优估计和相应的椭圆形置信 区域 其中已知的积分曲线满足常微分方程 0 0 0 dt ttT dt x m xxx 考虑到积分曲线的简单性及可解性 我们用由三角函数构成的且函数值在 0 和 1 之间的 简单函数作为已知的积分曲线 2 1 sin 8 3 2 1 cos 8 3 2 1 ttx ttx 原问题的已知条件为 12 1 ijijijijij YYi jn Ym xe 其中是的独立同分布的有界随机向量 为中的 12 1 n ijijij i j ee e 0 ij E e 12 ijijij xx x 2 1 0 网格点即 12 1 1 ijijij xxinjn x 我们假设已知 时 1 1 21 n j tx n i txnji 2 1 21 112 12 12 212 0 5 0 5 T T ij ij x tem x tx t x tx t x tem x tx t m 且有常微分方程 即 2 1 8 7 0 0 xxxm x t dt td 1 2 2 1 012 0 5 0 5 7 1 0 0 8 2 d x t x t dt d x t x t dt xx x 那么函数即是满足上述常微分方程的一个解 2 1 sin 8 3 2 1 cos 8 3 2 1 ttx ttx 我们利用MATLAB中的mvnrnd函数 Mvnrnd mu sigma cases 产生一个 cases d的服从多元正态分布的随机矩阵 其中均值向量mu是 1 d 协方差矩阵sigma是 d d的 但由于mvnrnd函数产生的值的随机性 我们事先取定一个随机种子数 以保证所 得结果可以被重复验证 已知的数据部分 程序如下 nji ijij 1 Yx randseed 100 randn seed randseed 苏州大学本科生毕业设计 论文 10 epi mvnrnd 0 0 0 1 0 0 0 1 n 2 for j 1 n for i 1 n x1 i n 1 x2 j n 1 Y i n j 1 1 x2 0 5 epi i n j 1 1 Y i n j 1 2 x1 0 5 epi i n j 1 2 end end 真实的积分曲线的值可直接由曲线的表达式求出 x t 0 pi 660 pi xx 3 8 cos x t 0 5 yy 3 8 sin x t 0 5 所以 主要问题在于对的求解 tx 4 24 2 积分曲线的最优估计及椭圆置信区域的求解积分曲线的最优估计及椭圆置信区域的求解 4 2 14 2 1 B 样条基函数矩阵样条基函数矩阵m文件的建立文件的建立 按照 3 1 中的B 样条函数的定义 若我们已知结点向量为 则根据B 样条函数的递归t 性质可求出第个m次 B 样条函数在向量上各点处的值 从而可以建立B 样条基函数矩jx 阵 该矩阵大小为numel numel n n代表样条为n 1次的 每一列为一个基函xt 数 函数名为bspline basismatrix n t x 有关M文件的代码见附录 4 2 24 2 2 的求解的求解 x m 根据 4 1 已知条件为 时 1 1 21 n j tx n i txnji 2 1 21 112 12 12 212 0 5 0 5 T T ij ij x tem x tx t x tx t x tem x tx t m 其中是的独立同分布的有界随机向量 我们采用三次B 样条的方 12 1 n ijijij i j ee e 0 ij E e 法来求解 x m 苏州大学本科生毕业设计 论文 11 首先用 4 2 1 中的m文件产生 B 样条基函数矩阵 函数bspline basismatrix的三XB 个自变量分别为p 4 结点向量 向量 11111 1 1 1 10000 NNNNXknots 其中n即为已知条件中的n N的选取满足1 1 1 1 1 11 nnnXx 此文中我们取 由此 我们得到了 2 1 4 4 221 5 p nNn log log 1 2 2 floornpnN 中所说的空间 这里 的一组基 mS1 pmXknots xB mj 在向量各分点处的值 B 样条基Nppj 2 1 1 1 1 1 11 nnnnXx 函数矩阵的具体矩阵表示为 XB 1 2 1 2 111 111 111 pmpmN m pmpmN m nNp BBB nnn nnn BBB nnn 定义张量积空间 1 2 21 1 0 n mjj mjjmjj xxBbmSmSxx 其中 张量积样条空间的维数为 2 12 1212 0 1 j m j j mj j mj m BBx xBx Bxx x x 由空间的定义 显然 基即为 2 pN 12 12 j m j j mj j mj m BBx xBx Bx xNppjj 2 1 对于基函数矩阵 用卡式积函数kron可求得双变量的 B 样条基函数矩阵 大小为XBA 即我们知道了每个双变量基函数在 2 pNn 12 j k mj k m BBx x x 处的取值 双变量的 B 样条基函数矩阵的具体矩 1 1 21 n k n j xxxnkj 2 1 A 阵表示形式为 1 1 1 2 1 1 1 2 1 1 1 2 111111 111111 121212 111111 111 pmpmpmpmN mN m pmpmpmpmN mN m pmpmpmpm BBBBBB nnnnnn BBBBBB nnnnnn nnnn BBBB nnnn 22 111 N mN m nNp nn BB nn 与一维的情况类似 其中的均可由张量积空间的 112 12 212 T m x x x x m x x m xm 21 m m 基线性表示 即 1 21 1 211 n mjj mjjjj xxBaxxm 1 21 2 212 n mjj mjjjj xxBaxxm 这样 求的问题 就归结为求上两式中的系数问题 求系数一般归结为解矩阵方程 m x 苏州大学本科生毕业设计 论文 12 4 1 节中我们已经产生了一组数据 其中 nji ijij 1 Yx 1 1 21 n j n i xx ij x 记 则求解系数归结为求解矩阵方程 12 ijijij YY Y 2 12 2 ijij n YY Y Y AB 其中 是随机误差矩阵 为已经求得的双变量 B 样条基函数矩阵 A 2 12 2 j jj j Np aa B 即为所求的系数矩阵 根据应用回归分析中多元线性回归的理论 我们用普通的最小二乘估 计来计算出回归系数矩阵 其中最小二乘法的过程详见参考文献 5 我们得到正规方程B 组 移项得 当存在时 即得回归系数矩阵的最小二乘估计 A YAB0 A ABA Y A A 为 则得到的样条估计在处的估计 BA AA Y m x x m 1 1 21 n j n i xx ij x 2 12 12 12 12 2 1111 1111 11 1111 2121 1111 1111 n mm nnnn nn mm nnnn mm nnnn nnnn mm nnnn 于是 我们得到向量场的样条估计的一般形式 写成矩阵形式即为 m x x m 为双变量 B 样条基函数在处取值产生的 12 TTT x x m xBA AA Y 21 xx T B 21 xx x 矩阵 这部分的代码为 n 140 p 4 N floor 1 n 2 2 p 1 log log n Xknots 0 0 0 0 1 N 1 1 N 1 N N 1 1 1 1 1 XB2test zeros 2 N p A zeros n 2 N p 2 XB bspline basismatrix p Xknots X1 A zeros n 2 N p 2 苏州大学本科生毕业设计 论文 13 for j 1 n for k 1 n A j 1 n k kron XB k XB j end end Vpin inv n 2 A A Beta n 2 Vpin A mhatbeta Beta Y mhatX A mhatbeta 4 2 34 2 3 的求解的求解 tx 我们采用欧拉方法来解微分方程 用已经求得的 0 0 0 dt ttT dt x m xxx m x 估计代替 已知 用公式 x m 2 1 8 7 0 0 xx 11 1 kkk tttkT xxm x 来一步步迭代 这里我们取步长 用 4 2 2 中求得的估计式05 0 k tm x 12 TTT kkk tx tx t m xBA AA Y 来计算 注意此时 为双变量B 样条基函数在处取值产生 21kk T txtxB 21kk txtx x 的矩阵 即需要再次调用 4 1 节中的m文件 这部分的代码如下 delta 0 005 ntest 628 xtest zeros 2 ntest xtest 1 7 8 1 2 XB1test zeros 2 N p ntest XB 1 zeros 2 N p 1 ntest XknotsP 0 0 0 1 N 1 1 N 1 N N 1 1 1 1 for k 1 ntest 1 if xtest 1 k 0 else XB1test 2 k zeros N p 1 end B0 kron XB1test 1 k XB1test 2 k mhatest B0 mhatbeta xtest k 1 xtest k B0 mhatbeta delta end 4 3 协方差矩阵以及椭圆置信区域的求解协方差矩阵以及椭圆置信区域的求解 4 3 14 3 1 协方差矩阵协方差矩阵的求解的求解 tG 对维矩阵的理论部分 我们这里只是做简单的公式介绍 由于 2 2 1 2 200 t t G 此部分内容超出本科生的知识范围 我们只是根据公式做相应的编程 对于推导跟证明过程 这里省略 的递推公式为 1 1 2 200 199 kkk k t ttk G 11111 kkkkpkkk t ttttptt GGmXXX 111111 T pkkkkkpk tttttt mXGGmX 其中 0 001 4p 1 12 TTT pkkk tx tx t mxBA AA Y 21 121 22 T ppp h X XBX VBX 其中 2 1 2 N p j j pjjpj jjjp n BB V 苏州大学本科生毕业设计 论文 15 其中 2 2 1 1111 n n i j ijij n nnnn 此部分的程序较长 我们直接附在附录里面 4 3 24 3 2 椭圆置信区域的求解椭圆置信区域的求解 根据定理 即 1 2 0 nhttNt XXG 11 22 10 0 01 nhtttN GXX 则根据正态分布和分布的关系有 2 1111 2 2222 2 TT nhtttnhttt XXGGXX 取 则的置信度为的椭圆置信区域的求解如下 05 0 tX 1 212 0 95 2 0 95 T P n httttt XXGXX 令 其中 12 ta a X 12 tx x X 1 1112 2 2122 bb t bb G 11 aa t 则 22 aa t 11 xx t 22 xx t 11 22 ax tt ax XX 由化简得 212 0 95 2 T n httttt XXGXX 2 1111122222 0 95 21112222 2 baxbax n h baxbax 2222 11111222211122220 95 2 baxbaxbaxbaxn h 这是一个椭圆的一般方程形式 令 即椭圆参数方程 2 0 95 11111222 2 2 0 95 21112222 2 2 cos 2 sin baxbax n h baxbax n h 02 用矩阵形式即为 2 0 95 1 2 2 2 0 95 2 2 cos 2 sin n h ttt n h GXX 则 2 0 95 2 11 2 0 95 22 1 22 2 0 950 95 2 2 cos 2 cos 1 2 2 sin sin n h ttttt nh n h XXGXG 程序如下 苏州大学本科生毕业设计 论文 16 alpha 0 05 r sqrt chi2inv 1 alpha 2 theta 0 0 05 1 2 pi len length theta z zeros len 2 2 for j 1 len can change to another method z j 1 r cos theta j z j 2 r sin theta j end for k 1 ntest 1 if abs k 1 15 floor k 1 15 10 eps ellipse repmat xtest k 1 1 len real sqrtm var est k 1 z n N 1 0 5 plot3 ellipse 1 ellipse 2 k 1 ones 1 len b LineWidth 1 其中 我们利用循环语句画出置信椭圆区域 但是如果每一个k处都画一个椭圆 椭圆太密 在图上虽然能细致的看出置信椭圆带的轮廓 但椭圆带内的估计的积分曲线和真实的积分曲 线根本看不清 如下图 所以我们考虑只选取660个椭圆中的一部分椭圆 因此增加了条件 苏州大学本科生毕业设计 论文 17 if abs k 1 15 floor k 1 15 10 eps 其中eps是一个非常小的数 只有满足此 条件的k 我们才会在图上画椭圆 画的椭圆个数是660 15 44个 例如上图重新作图之后 这样的话 图上的椭圆以及已知积分曲线和它的估计曲线都比较清晰 而且也不会影响整个 椭圆置信带的布局结构 第五章第五章 结果分析结果分析 本文的目的是得到已知积分曲线的最优估计 以及置信椭圆带能较大程度包括住积分曲 线 所以在已知已知积分曲线时 我们可根据实际的积分曲线来判断估计的好坏 5 5 1 1 对积分曲线整体估计好坏的简单分析对积分曲线整体估计好坏的简单分析 首先不同的n 不仅计算量不同 而且置信椭圆带的好坏 即是否直观上很好的把真实曲 线包括在里面 也会不同 下面我们就根据选取的已知积分曲线来改进程序 上面我们已求出估计值和置信椭 t x 圆 且均画出了置信椭圆带和估计的积分曲线 一方面我们可以直接观察置信椭圆带的宽度 另外 可对每个k 判断已知积分曲线的真实值是否在相应的置信椭圆内 看置信椭圆带包 苏州大学本科生毕业设计 论文 18 括住真实的积分曲线的程度大小 我们可得出660组数据 有多少组真实值在相应置信椭圆 内 这部分代码如下 n0 0 for k 2 660 if inv sqrtm var est k xtest k xx k yy k inv sqrtm var est k xtest k xx k yy k 250 越来越大 660 组真实值在置信椭圆中的数目 越来越少 综合考虑结果和计算速度 我们得出 n 在 100 至 150 之间 置信椭圆带能够较大 程度的包括住已知积分曲线 且计算速度较快 5 5 2 2 基于置信度概念的估计好坏的判断基于置信度概念的估计好坏的判断 在 4 1 节中我们取定种子数 用mvnrnd函数来产生服从多元正态分布的随机矩阵 程 序中我们只取了一组随机数 也就是的值是只取了一组 12 1 ijijijijij YYi jn Ym xe 样本 得到了的一个估计以及相应的置信椭圆带 但是对于一次抽样 判断得到的椭圆 tX 置信区域是否包括住在处的真实值是没有统计意义的 基于置信度的概念 我们需要 tX k t 重复取样 对固定的一个 会得到许多不同的椭圆置信区域 我们看实际上这些椭圆置信 k t 区域中包括住在处的真实值有多少 如果能够达到 95 在 4 3 节中我们取的置信度 tX k t 为 95 则说明我们的估计方法比较好 首先我们为了取样 需要产生 100 组 即需要随机 12 1 ijijijijij YYi jn Ym xe 数有 100 组 这部分程序如下 苏州大学本科生毕业设计 论文 20 n 100 randseed 100 randn seed randseed x mvnrnd 0 0 0 1 0 0 0 1 n 2 100 epi zeros n 2 2 100 epi 1 x 1 10000 for i 2 100 epi i x i 1 10000 1 i 10000 end for a 1 100 for j 1 n for i 1 n x1 i n 1 x2 j n 1 Y i n j 1 1 a x2 0 5 epi i n j 1 1 a Y i n j 1 2 a x1 0 5 epi i n j 1 2 a end end 至此 我们已经产生了 100 组样本 我们对每组样本都进行 4 1 4 4 节已经编好的程序计 算 得到每组样本所对应的 660 个椭圆形置信区域 对每个 我们判断 100 个样本得到的 k t 置信椭圆 有多少个能包括住在处的真实值 tX k t 这部分 我们把原来程序中的增加一维即程序中xtest为一个三维的矩阵 第三维 t X 代表样本 在对k的循环里 加入对样本的循环 此时 我们在循环里判断 100 个样本得到 的置信椭圆是否包括住在处的真实值 包括住的数目为n1 得出 100 个样本中包括住 tX k t 真实值的比率 矩阵b即为 660 个n1 的值 所以这就使得程序的计算量为原来的 100 倍 一部分改动的程序如下 xtest zeros 2 ntest 100 for a 1 100 xtest 1 a 7 8 1 2 end b zeros 660 1 for k 1 ntest 1 n1 0 苏州大学本科生毕业设计 论文 21 for a 1 100 这部分程序与对一个样本做估计时的程序相同 这里省略 详细见附录 if inv sqrtm var est k 1 xtest k a xx k yy k inv sqrtm var est k 1 xtest k a xx k yy k 1 b bspline basis j n 1 t x dn x t j 1 dd t j n t j 1 if dd 0 indeterminate forms 0 0 are deemed to be zero y y b dn dd end b bspline basis j 1 n 1 t x dn t j n 1 x dd t j n 1 t j 1 1 if dd 0 y y b dn dd end elseif t j 2 t end treat last element of knot vector as a special case y t j 1 x else y t j 1 2 B zeros numel x numel t n for j 0 numel t n 1 B j 1 bspline basis j n t x end else b x bspline basis 0 n t B zeros numel x numel t n B 1 b for j 1 numel t n 1 B j 1 bspline basis j n t x end end 已知已知积分曲线 积分曲线 randseed 100 randn seed randseed n 140 p 4 N floor 1 n 2 2 p 1 log log n 苏州大学本科生毕业设计 论文 25 Xknots 0 0 0 0 1 N 1 1 N 1 N N 1 1 1 1 1 XB2test zeros 2 N p A zeros n 2 N p 2 epi mvnrnd 0 0 0 1 0 0 0 1 n 2 for j 1 n for i 1 n x1 i n 1 x2 j n 1 Y i n j 1 1 x2 0 5 epi i n j 1 1 Y i n j 1 2 x1 0 5 epi i n j 1 2 end end ellipse plot set up alpha 0 05 r sqrt chi2inv 1 alpha 2 theta 0 0 05 1 2 pi len length theta z zeros len 2 for j 1 len z j 1 r cos theta j z j 2 r sin theta j end estimate Sigma matrix X1 1 n 1 1 n 1 n n 1 XB bspline basismatrix p Xknots X1 spline basis matrices A zeros n 2 N p 2 for j 1 n for k 1 n A j 1 n k kron XB k XB j end end Vpin inv n 2 A A Beta n 2 Vpin A mhatbeta Beta Y mhatX A mhatbeta Sigma hat Y mhatX Y mhatX n 2 generate new Vp for estimating rho function new 200 Nnew floor new 2 2 p 1 log log new Xknotsnew 0 0 0 0 1 Nnew 1 1 Nnew 1 Nnew Nnew 1 1 1 1 1 X1new 1 new 1 1 new 1 new new 1 XBnew bspline basismatrix p Xknotsnew X1new Anew zeros new 2 Nnew p 2 苏州大学本科生毕业设计 论文 26 for j 1 new for k 1 new Anew j 1 new k kron XBnew k XBnew j end end Vpinnew inv new 2 Anew Anew initial steps for estimate the derivative of m and x t delta 0 005 ntest 660 xtest zeros 2 ntest xtest 1 7 8 1 2 D eye N p 1 diag ones 1 N p 1 1 A1 D 1 N p 1 s 1 for i 1 p s 1 A1 i A1 i p s p p s i 1 A1 N p i s 1 A1 N p i s 1 p s p p s i 1 end D1 A1 XB1test zeros 2 N p ntest XB1test new zeros 2 Nnew p ntest XB 1 zeros 2 N p 1 ntest mhat v zeros 2 2 ntest var est zeros 2 2 ntest