• 沒有找到結果。

在僅能觀測到部分資料的情形下,ODE model的參數估計

N/A
N/A
Protected

Academic year: 2021

Share "在僅能觀測到部分資料的情形下,ODE model的參數估計"

Copied!
35
0
0

加載中.... (立即查看全文)

全文

(1)國立臺灣師範大學數學系碩士班碩士論文. 指導教授:陳賢修 博士. Parameter Estimation for Partially Observed ODE and its Applications. 研 究 生:彭瑋翔. 中 華 民 國 一百零二 年 二 月.

(2) 感謝詞 首先我要感謝我的指導教授. 陳賢修老師,由於他兩年. 來耐心的指導,以及在我遇到瓶頸的時候,給予我許多有建 設性的建議和幫助,因此我才能順利完成這篇論文。在此由 衷感謝老師。 此外我還要感謝. 鄭昌源老師 和 黃聰明老師百忙之. 中抽空擔任口試委員,並且對文章提出建議,指出不足之缺 陷,使得論文能夠更加完備。 最後要感謝研究室的各位同學,這兩年一同陪伴、一起 歡笑的日子,讓我做論文之餘也能感受到你們的關心與鼓 勵。 中華民國 102 年 2 月.

(3) 中文摘要. 在物理、化學、生物的領域上,我們常常會用微分方程 模型去描述自然的現象。而對於這樣一個微分方程模型,我 們可能僅能觀測到裡面的其中某幾個參數。在這樣的情形底 下,我們該如何去做微分方程模型的參數估計?這就是我們 要討論的議題。在考慮到電腦計算的速度以及參數估計的精 確度,本文以 Hindmarsh-Rose Type model 做為例子,提供了 一套參數估計的方法。 關 鍵 字:參數估計、微分方程模型、類 Hindmarsh-Rose 模 型.

(4) Parameter Estimation for Partially Observed ODE and its Applications. February 26, 2013. Abstract In physics, chemistry and biology, many researcher try to describe natural phenomena by ordinary differential equations (ODE) due to the fact that a lot of elegant theorem can be found in ODE and its nonlinear properties. With the observed information and some prior knowledge for a potential model, parameter estimation of an ODE model often requires a lot of computation works, such as a numerical integration of the ODE system and minimization of the log-likelihood function. In this paper, we consider a partially observed ODE model: 2D-HR type model. With the partially information, we present a simple way (iterative estimation process) without too mach computation to estimate the parameters in this model.. keywords: parameter estimation, ordinary differential equations, Hindmarsh-Rose Type model, partially observed model.. 1.

(5) Contents 1 Introduction. 3. 2 Hindmarsh-Rose model and our main problem. 4. 3 Ideas for Parameter Estimation on 2D-HR type model. 5. 4 simulation on 2D-HR type model without noise. 8. 5 simulation on 2D-HR type model with Gaussian noise. 20. 6 Conclusion. 25. 7 Some Tables 7.1 The simulation result while ∆t decreases with k=40, the initial c0 = 10 7.2 The simulation result while ∆t decreases with k=80, the initial c0 = 10 7.3 The simulation result while ∆t decreases with k=40, the initial c0 = 3 .. 26 26 28 30. 8 Reference. 32. 2.

(6) 1. Introduction. In physics, chemistry and biology, ordinary differential equations (ODE) is a wide used model to describe natural phenomena. For a ODE model, we also need to understand the feature by estimation of parameters from the information and observations. Therefore, parameter estimation is an important issue in many fields of science. In recent years, Scientists have developed many method of parameter estimation on a nonpartially observed ODE model. Usually, they may estimate the parameters by classical parametric estimators, like the least squares estimator(brunel[1]) and Maximum likelihood estimator(MLE), or by other different method. However, despite their satisfactory theoretical properties, the efficiency of these estimators may be dramatically degraded in practice by computational problems that arise from the implicit and nonlinear definition of the model. Indeed, these estimators give rise to nonlinear optimization problems that necessitate the approximation of the solution of the ODE and the exploration of the (usually high-dimensional) parameter space. Hence, we have to face possibly numerous local optima and a huge computation time. Under these problem, we need to find a method of simple computations which can estimate the parameters of a partially observed ODE. In this paper, we develop a iterative estimation process. We shows our iterative estimation process method and the simulation results by 2D-HR type model. In section 2, we introduce 2D-HR type model and our main problem. The idea of iterative estimation process is showed in section 3. By this process, We put the simulation results of parameter estimation in section 4,5. After all, we conclude the result in section 6.. 3.

(7) 2. Hindmarsh-Rose model and our main problem. In neuron science, Hindmarsh-Rose model(HR-model) is composed by ordinary differential equations. This model is aimed to study and to describe the behavior of the membrane potential observed in experiments with a single neuron. In order to study the parameter estimation problem on partially observed ODE, in this paper, we will focus on a simplified HR-model: 2-dimension Hindmarsh-Rose type model(2DHR type model): x3 − y), 3  x2 + dx − by + a /c.. x˙ = c(x −. (2.1). y˙ =. (2.2). This model contains two variables: the membrane potential of a single neuron x and the recovery variable y. In experiments, we can measure the membrane potential x of neurons by machines, and get a discrete data of x with time. But, the recovery variable y is a auxiliary variable to describe the membrane potential x, and can not be observed by machines. Here, we can not observed the data all of the variables, so this 2D-HR type model is a partially observed model. In the following, our main problem is how to estimate the parameters of 2D-HR type model. Then, we want to extend this results to parameter-estimation on other partially observed model.. Figure 2.1: Simulation data of the membrane potential: simulated by 2D-HR type model where a = 1.05, b = 2, c = 3, and d = 3.55. X-axis is the time t and Y-axis is the membrane potential x. 4.

(8) 3. Ideas for Parameter Estimation on 2D-HR type model. How to estimate the parameters(a, b, c, d) in 2D-HR type model where the data of x(t) is observed and y(t) is not? We will show some ideas in this section. In 2D-HR type model ( (2.1), (2.2)), our first trouble is that y(t) can not be observed. To solve this problem, we write (2.1) to another form: y = x−. x3 x˙ − , 3 c. (3.1). and give an initial value of c (we call this initial c0 ), then we can calculate the value of y(t) with time by x(t): y0 = x −. x3 x˙ − . 3 c0. (3.2). Help by (3.2), y(t) can be calculated by x(t) and c0 , and we named it y0 (t). Because x is discrete, our estimation results will depend on the choose of calculation of dotx(t). We discuss this problem in section 5.. (a) c0 = 3. (b) c0 = 10. Figure 3.2: Figure of y0 (y-axis) with time (x-axis): Calculated by different c0 and the simulation data of x in fig 2.1. 5.

(9) Because this y0 (t) and c0 is a temporary answer for y(t) and c, next we want to modify y and c to approach the real value. Rewrite (2.2) as following : y˙ =. 1 2 d b a x + x− y+ , c c c c. (3.3). and we regard this equation to a normal form: Y = Xβ. (3.4). where Y = (y(t))T , X = (x2 (t), x(t), y(t), 1)T , and the parameter of interest is β = (β1 , β2 , β3 , β4 ) 1 d b a = ( , , − , )T . c c c c. (3.5) (3.6). We have a least-squares estimator: βˆ = (X T X)−1 X T Y.. (3.7). ˆ We also get an estimator for a, b, c, and d: From β, β4 β1 ˆb = − β3 β1 1 cˆ = β1 β2 dˆ = β1. a ˆ =. (3.8) (3.9) (3.10) (3.11). With y0 we computed by c0 , and the observation of x, we can estimate the parameters in our 2D-HR type model by the equations above, and we named the value we estimate here by a1 , b1 , c1 , and d1 , respectively. With this new value of c: c1 , and the observation of x( t), we can use (3.1) to get a new estimate of y(t)(named y1 (t)). Then, repeat this process, we can also get a new 6.

(10) value c2 , c3 , c4 , and so on. With this iterative process, we expect that we can get a more exact estimator for the parameters (a, b, c, and d) in 2D-HR model. In the following, we sort out the steps of this iterative process: An iterative estimation process for partially observed 2D-HR type model: (Suppose the iterative times is k) Step1 Choose an initial value of c (c0 ) Step2 By (3.2), substitute ck and the observation of x(t), calculate yk (t). Step3 Use the observation of x, and yk (t) we calculate in previous step, estimate the parameters in 2D-HR model by (3.3). We get ak+1 , bk+1 , ck+1 , and dk+1 (least-squares estimator) Step4 Take the new value of c (ck+1 ), and repeat Step2.. Figure 3.3: The iterative process of parameter estimation for partially observed 2D-HR type model. In next section, we will show the accuracy of this iterative estimation process by simulation. 7.

(11) 4. simulation on 2D-HR type model without noise. In this section, we will show the estimation results by iterative estimation process with an simulation data. In partially observed 2D-HR type model, we suppose the solution of x(t) can be observed at time ti , i = 0, 1, ..., n, where 0 = t0 < ... < tn = T , and the time step ∆t = (ti+1 − ti ) = T /n is constant. Also, with this partially observed ODE, we suppose the solution of y is unknown. Now, to simulate the parameter of 2D-HR type model with this iterative process, we still need to compute the differentiation of x(t) (Because the observation of x(t) is discrete). Let ∆t be the time step of sampling, we compute dotx by the following equation: x˙ t+ 2t ∆t =. x(t + h) − x(t) , ∆t. (4.1). Also, we compute the value of x at time t + 2t h by the mean of x(t + h) and x(t): xt+ 2t ∆t =. x(t + ∆t) + x(t) . 2. (4.2). By the equtions (4.1) and (4.2), the value of x and x˙ at every time ti and ti+ 2t ∆t can be computed for all i = 0, 1, ..., n − 1, (This is a modification of central difference method. We compare the simulation results of our method and central difference in Appendix) From these data, we start to estimate the parameters of 2D-HR type model by our iterative estimation process, and show the simulation results as following. Accuracy of the estimator yˆk First, for the iterative estimation process, we need to estimate y by calculation on (3.2). In figure 4.5, we show the simulation results of the difference of the original solution y and yˆk while k = 1, 2, 3. Here we set a = 1.05, b = 2, c = 3, d = 3.55, a parameter setting near Hopf-bifurcation in Dynamic system. The initial value of c (c0 ) is 10, ∆t = 2−9 , and total time T = 50. We can see that yˆk (t) (the blue curve) close to y (the red curve) as k increased. This result tells us our iterative estimation process and can be well used to estimate y at this parameter setting.. 8.

(12) Accuracy of the estimators a ˆk , ˆbk , cˆk , and dˆk With a good estimate of y, we want to know the accuracy of the parameters at the parameter setting above. Table 4.1 show that the parameters we estimate by our iterative estimation process is also close to the true parameter as k increased. When k = 10, the difference of c and cˆk is 0.0187. When k = 20, the difference is only 0.0001. Furthermore, when k = 20, a ˆk , ˆbk , cˆk , and dˆk , also has difference 0.0001 with the true parameter. It shows that this iterative estimation process works good at this parameter setting. Furthermore, in the following, we also simulate at other parameter settings: Bautinbifurcation(a = 0.518, b = 1.5, c = 3, d = 2.6) and BT bifurcation(a = −0.09715, b = 2, c = 3, d = 2.2528). We show the result on the following table and diagram. Similarly, as iterative times k increased, we can estimate the parameter a, b, c, d accurately at these two parameter settings.. Figure 4.4: The bifurcations in 2D-HR type model while b = 2 and c = 3, The X-axis is a, and the Y-axis is d. The red curve is Hopf bifurcation, green and black curve are Saddle-node bifurcation. 9.

(13) (a) k = 1. (b) k = 2. (c) k = 3. Figure 4.5: Diagram of yˆk (t) (the blue curve) and y (the red curve) with k = 1, 2, 3 near Hopf-bifurcation, where a = 1.05, b = 2, c = 3, d = 3.55, T = 50, ∆t = 2−9 , and initial c0 = 10 The horizontal axis is time t, and the vertical axis is the value of y.. 10.

(14) Table 4.1: Simulation results while k increases: 2D-HR model with a = 1.05, b = 2, c = 3, d = 3.55, T = 50, ∆t = 2−9 , and initial c0 = 10. k 1 2 3 4 5 6 7 8 9 10. aˆk -0.7075 0.1486 0.5314 0.7416 0.8643 0.9375 0.9817 1.0084 1.0247 1.0346. bˆk 2.7519 2.3865 2.2247 2.1347 2.0815 2.0495 2.0302 2.0184 2.0112 2.0068. cˆk 5.5153 4.2497 3.6864 3.3948 3.2327 3.1390 3.0838 3.0507 3.0308 3.0187. dˆk 2.0883 2.8009 3.1200 3.2947 3.3964 3.4570 3.4936 3.5157 3.5291 3.5373. bias of cˆk and c 2.5153 1.2497 0.6864 0.3948 0.2327 0.1390 0.0838 0.0507 0.0308 0.0187. 11 12 13 14 15 16 17 18 19 20. 1.0406 1.0443 1.0465 1.0479 1.0487 1.0492 1.0495 1.0497 1.0498 1.0499. 2.0042 2.0025 2.0015 2.0009 2.0006 2.0004 2.0002 2.0001 2.0001 2.0001. 3.0114 3.0069 3.0042 3.0026 3.0016 3.0009 3.0006 3.0003 3.0002 3.0001. 3.5422 3.5453 3.5471 3.5482 3.5489 3.5493 3.5496 3.5497 3.5498 3.5499. 0.0114 0.0069 0.0042 0.0026 0.0016 0.0009 0.0006 0.0003 0.0002 0.0001. 11.

(15) (a) k = 1. (b) k = 2. (c) k = 3. Figure 4.6: Diagram of ykˆ(t) (the blue curve) and y (the red curve) with k = 1, 2, 3 near Bautin bifurcation, where a = 0.518, b = 1.5, c = 3, d = 2.6, T = 50, ∆t = 2−9 , and initial c0 = 10 The horizontal 12 axis is time t, and the vertical axis is the value of y..

(16) Table 4.2: Simulation results while k increases: 2D-HR model on Bautin bifurcation with a = 0.518, b = 1.5, c = 3, d = 2.6, T = 50, ∆t = 2−9 , and initial c0 = 10. k 1 2 3 4 5 6 7 8 9 10. aˆk -0.5735 -0.0576 0.2864 0.4383 0.4918 0.5095 0.5153 0.5171 0.5177 0.5179. bˆk 1.8294 1.7120 1.5904 1.5316 1.5104 1.5034 1.5011 1.5003 1.5001 1.5000. cˆk 4.9675 3.5891 3.1824 3.0579 3.0185 3.0059 3.0019 3.0006 3.0002 3.0000. dˆk 1.6707 2.1083 2.4019 2.5318 2.5776 2.5927 2.5977 2.5992 2.5997 2.5999. bias of cˆk and c 1.9675 0.5891 0.1824 0.0579 0.0185 0.0059 0.0019 0.0006 0.0002 0.0000. 11 12 13 14 15 16 17 18 19 20. 0.5180 0.5180 0.5180 0.5180 0.5180 0.5180 0.5180 0.5180 0.5180 0.5180. 1.5000 1.5000 1.5000 1.5000 1.5000 1.5000 1.5000 1.5000 1.5000 1.5000. 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000. 2.6000 2.6000 2.6000 2.6000 2.6000 2.6000 2.6000 2.6000 2.6000 2.6000. -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000. 13.

(17) (a) k = 1. (b) k = 2. (c) k = 3. Figure 4.7: Diagram of ykˆ(t) (the blue curve) and y (the red curve) with k = 1, 2, 3 near BT bifurcation, where a = −0.09715, b = 2, c = 3, d = 2.2528, T = 50, ∆t = 2−9 , and initial c0 = 10 The horizontal axis is time t, and the vertical axis is the 14 value of y..

(18) Table 4.3: Simulation results while k increases: 2D-HR model on BT bifurcation with a = −0.09715, b = 2, c = 3, d = 2.2528, T = 50, ∆t = 2−9 , and initial c0 = 10. k 1 2 3 4 5 6 7 8 9 10. aˆk -1.1159 -0.6015 -0.3300 -0.1977 -0.1392 -0.1145 -0.1043 -0.1001 -0.0983 -0.0976. bˆk 2.2820 2.1697 2.0849 2.0378 2.0160 2.0066 2.0027 2.0011 2.0005 2.0002. cˆk 4.9939 3.7244 3.2820 3.1129 3.0457 3.0186 3.0076 3.0031 3.0013 3.0005. dˆk 1.4055 1.8427 2.0655 2.1723 2.2192 2.2390 2.2471 2.2505 2.2518 2.2524. bias of cˆk and c 1.9939 0.7244 0.2820 0.1129 0.0457 0.0186 0.0076 0.0031 0.0013 0.0005. 11 12 13 14 15 16 17 18 19 20. -0.0974 -0.0972 -0.0972 -0.0972 -0.0972 -0.0972 -0.0972 -0.0972 -0.0972 -0.0972. 2.0001 2.0000 2.0000 2.0000 2.0000 2.0000 2.0000 2.0000 2.0000 2.0000. 3.0002 3.0001 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000 3.0000. 2.2526 2.2527 2.2528 2.2528 2.2528 2.2528 2.2528 2.2528 2.2528 2.2528. 0.0002 0.0001 0.0000 0.0000 0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000. 15.

(19) Other parameter settings. Now, we want to know how the iterative estimation process works at other parameter settings. In figure 4.8, we fixed b = 2, c = 3, and let a changes from -2 to 1, d changes from 1 to 3. Also, let the iterative times k = 40, T = 50, and ∆t = 2−9 . In this figure, The bias of cˆk and c in blue area is less than 0.001. in green area, the bias is between 0.1 and 0.001. In red area, the bias is greater than 0.1.. Figure 4.8: Simulation results: parameter estimation by iterative estimation process when b = 2 and c = 3. X-axis is the true value of a and Y-axis is the true value of d.. 16.

(20) (a) ck and ck+1 at a = 0 d = 2.5 (b) dcdk ck+1 at a = 0 d = 2.5 in in blue area blue area. (c) ck and ck+1 at a = −1.1 d = (d) dcdk ck+1 at a = −1.1 d = 2.4 2.4 in green area in green area. (e) ck and ck+1 at a = −1.3 d = (f) dcdk ck+1 at a = −1.3 d = 2.6 2.6 in red area in red area. Figure 4.9: Diagram of ck and ck+1 and in blue, green, and red area: fixed b = 2, c = 3, T = 50, and dt = 2−9 left-hand side: the X-axis is ck , the Y-axis is ck+1 . right-hand side: the X-axis is ck , the Y-axis is dcdk ck+1 .. 17.

(21) Here, we found that there is some region that iterative estimation process can not be well used. In order to find the reasons of the difference of red, green, and blue area, we plot figure 4.9 to see how the iterative estimation process works. In the lefthand side of figure 4.9, X-axis is the value of ck , and Y-axis is the value of ck+1 . In the right-hand side, X-axis is the value of ck , and Y-axis is the differential of ck+1 ( dcdk ck+1 ). By Newton iterative method, if the differential of ck+1 is less then 1 in some neighborhood containing the true parameter, then our iterative estimation process for 2D-HR type model will converges in this neighborhood. From figure 4.9, we could see that in the blue area, dcdk ck+1 is smaller then 1 near c = 3. so parameter ck could converges to the true value in this area. In the green area, c˙k is also smaller then 1 near c = 3, but compare with the blue area, dcdk ck+1 is relatively close to 1 near c = 3. This means that ck in the green area will converge to the true value slowly. In the red area, because d c is greater then 1 at c = 3, the value of c will diverge. dck k+1 Different time steps ∆t. (figure 4.10) In the following, we show that how the time step ∆t affect our simulation results. With different time step ∆t (2−4 , 2−5 , ..., 2−9 ), we consider −2 ≤ a ≤ 1, b = 2, c = 3, and1 ≤ a ≤ 3. In figure 4.10, the area of blue region increased as ∆t decreased. This result is quite intuitive. As the number of observed data increased, the accuracy of parameter estimation will also increased. Similarly, the area of green region is also increased when ∆t decreased. This shows that the effect of estimation is better when ∆t is small.. After these simulation work, we could see that there’s three factor can effect our estimation results: the iteration times k, the initial value c0 , and the size of ∆t. In Section 7, we show some tables with different k, c0 , and ∆t. From these table, we know that when | c0 − c | and ∆t are sufficiently small, our iterative estimation process can be well used.. 18.

(22) (a) ∆t = 2−4 ≈ 0.0625. (b) ∆t = 2−5. (c) ∆t = 2−6 ≈ 0.016. (d) ∆t = 2−7. (e) ∆t = 2−8. (f) ∆t = 2−9 ≈ 0.002. Figure 4.10: Diagram of diversity for different time steps with the fixed parameter b = 2, c = 3. As dt decreased, the accuracy will increase, too.. 19.

(23) 5. simulation on 2D-HR type model with Gaussian noise. In this section, we want to see the availability of this iterative estimation process when the observations have a Gaussian noise: Xobs (ti ) = x(ti ) + i , i = 0, 1, ..., n. (5.1). where x(ti ) is the solution of x in 2D-HR type model at time t, and i follows normal distribution independently for each i. In order to increase the accuracy of estimation, we need to reduce the noise from the observation. Here we choose a smooth method: Nadaraya-Watson estimator, which is a nonparametric kernel regression estimator and can smooth the observation. Now, we estimate the value of x at each time t on the basis of Nadaraya-Watson estimator. Let Pn t−ti i=0 K( h )Xobs (ti ) Pn , (5.2) xˆ(t) = t−ti i=0 K( h ) where h is the bandwidth, and K is a kernel function with Gaussian kernel: 1 2 K(z) = √ e−z . 2π. (5.3). To use Nadaraya-Watson estimator, we need to choose a good bandwidth h before smoothing. In table 5.4 and 5.5, we set the parameters in 2D-HR type model near BT bifurcation: a = −0.09715, b = 2, c = 3, d = 2.2528, dt = 2−7 , and T = 50. These two table show the maximum difference between xˆ(ti ) and x(ti ), and the maximum difference between xˆ˙ (ti ) and x(t ˙ i ). The white noise level 10% means the energy i of noise is the tenth of observation, 20% means the energy i is the fifth of observation, and so on. In this table, each result(mean and standard error) of different white noise level and bandwidth h is calculated by repeat 50 times.. 20.

(24) Table 5.4: Find the best h in Nadaraya-Watson estimator by max | x ˆ(ti ) − x(ti ) |: 2D-HR model with different white noise level repeat 50 times, where a = −0.09715, b = 2, c = 3, d = 2.2528, ∆t = 2−7 , and T = 50. 5% white noise max | xˆ(ti ) − x(ti ) |. max | xˆ(ti ) − x(ti ) |. h 0.001 0.002 0.005 0.01 0.02 0.05 0.1 0.2 0.5 1 2 5. h 0.001 0.002 0.005 0.01 0.02 0.05 0.1 0.2 0.5 1 2 5. Mean 0.0335 0.0335 0.0332 0.0460 0.0838 0.1939 0.3432 0.5229 0.5969 1.0057 1.7088 2.5534. SD 0.0000 0.0001 0.0015 0.0053 0.0045 0.0029 0.0020 0.0015 0.0009 0.0004 0.0004 0.0003. 10% white noise. 20% white noise. Mean 0.0670 0.0670 0.0656 0.0664 0.1003 0.2096 0.3605 0.5396 0.6138 0.9888 1.6922 2.5366. Mean 0.1341 0.1340 0.1301 0.1217 0.1315 0.2446 0.3928 0.5728 0.6475 1.0070 1.6587 2.5032. SD 0.0000 0.0000 0.0017 0.0074 0.0087 0.0053 0.0038 0.0026 0.0016 0.0010 0.0007 0.0006. SD 0.0000 0.0000 0.0015 0.0047 0.0153 0.0101 0.0081 0.0055 0.0034 0.0015 0.0012 0.0009. 30% white noise. 40% white noise. 50% white noise. Mean 0.2011 0.2010 0.1958 0.1836 0.1798 0.2700 0.4260 0.6062 0.6816 1.0405 1.6251 2.4697. Mean 0.2681 0.2680 0.2605 0.2433 0.2241 0.3077 0.4620 0.6393 0.7137 1.0728 1.5922 2.4359. Mean 0.3352 0.3350 0.3249 0.3028 0.2752 0.3416 0.4949 0.6763 0.7488 1.1080 1.5589 2.4019. SD 0.0000 0.0000 0.0022 0.0076 0.0177 0.0190 0.0130 0.0074 0.0059 0.0028 0.0016 0.0016. 21. SD 0.0000 0.0001 0.0025 0.0056 0.0184 0.0245 0.0174 0.0124 0.0066 0.0035 0.0025 0.0023. SD 0.0001 0.0001 0.0030 0.0083 0.0127 0.0281 0.0215 0.0137 0.0083 0.0040 0.0031 0.0031.

(25) ˆ˙ i ) − x(t Table 5.5: Find the best h in Nadaraya-Watson estimator by max | x(t ˙ i ) |: 2D-HR model with different white noise level repeat 50 times, where a = −0.09715, b = 2, c = 3, d = 2.2528, ∆t = 2−7 , and T = 50. 5% white noise ˆ˙ i ) − x(t max | x(t ˙ i) |. ˆ˙ i ) − x(t max | x(t ˙ i) |. h 0.001 0.002 0.005 0.01 0.02 0.05 0.1 0.2 0.5 1 2 5. h 0.001 0.002 0.005 0.01 0.02 0.05 0.1 0.2 0.5 1 2 5. Mean 0.0331 0.0330 0.0194 0.0175 0.0218 0.0252 0.0286 0.0336 0.0392 0.0420 0.0419 0.0396. SD 0.0002 0.0002 0.0012 0.0023 0.0009 0.0002 0.0001 0.0000 0.0000 0.0000 0.0000 0.0000. 10% white noise. 20% white noise. Mean 0.0664 0.0662 0.0381 0.0191 0.0219 0.0252 0.0286 0.0336 0.0392 0.0420 0.0419 0.0396. Mean 0.1326 0.1323 0.0757 0.0324 0.0223 0.0250 0.0285 0.0336 0.0392 0.0420 0.0419 0.0396. SD 0.0004 0.0004 0.0010 0.0038 0.0019 0.0005 0.0002 0.0001 0.0000 0.0000 0.0000 0.0000. SD 0.0008 0.0009 0.0017 0.0017 0.0038 0.0008 0.0004 0.0001 0.0000 0.0000 0.0000 0.0000. 30% white noise. 40% white noise. 50% white noise. Mean 0.1988 0.1988 0.1141 0.0482 0.0242 0.0251 0.0286 0.0336 0.0392 0.0421 0.0419 0.0396. Mean 0.2646 0.2647 0.1517 0.0642 0.0269 0.0250 0.0287 0.0336 0.0392 0.0420 0.0419 0.0396. Mean 0.3312 0.3308 0.1900 0.0808 0.0316 0.0252 0.0286 0.0336 0.0392 0.0420 0.0419 0.0396. SD 0.0013 0.0011 0.0034 0.0020 0.0046 0.0011 0.0004 0.0001 0.0000 0.0000 0.0000 0.0000. 22. SD 0.0017 0.0016 0.0043 0.0030 0.0037 0.0018 0.0007 0.0002 0.0001 0.0000 0.0000 0.0000. SD 0.0019 0.0017 0.0051 0.0040 0.0029 0.0022 0.0007 0.0003 0.0001 0.0000 0.0000 0.0000.

(26) In this table, we consider the value of h in different white noise i such that max | xˆ(ti ) − x(ti ) | and max | xˆ(ti ) − x(ti ) | are small. The variance of this two max difference also can not be too large. Combine our test result, we choose h = 0.02 when white noise level is 5%, h = 0.05 when white noise level is 10%, 20%, 30%. h = 0.1 when white noise level is 40%. parameter estimation from noisy observation. Here, we choose the parameter setting near Bautin bifurcation, a = −0.518, b = 1.5, c = 3, d = 2.6, T = 50, ∆t = 2−7 , in different white noise level. Table 5.6: Simulation results while white noise level increases: 2D-HR model on BT bifurcation with a = −0.518, b = 1.5, c = 3, d = 2.6, T = 50, ∆t = 2−7 , and initial c0 = 10. white noise level 0% h a ˆ ˆb cˆ dˆ. Mean SD Mean SD Mean SD Mean SD. 0.5179 1.4998 2.9992 2.5999. 5% 0.01 0.5148 0 1.4910 0 2.9814 0 2.5964 0. 10% 0.02 0.4754 0.0047 1.6295 0.0156 3.1186 0.0275 2.6268 0.0045. 20% 0.05 0.3898 0.0028 1.5925 0.0116 2.8468 0.0210 2.5668 0.0040. 30% 0.05 0.2812 0.0066 1.8410 0.0305 3.0283 0.0550 2.6025 0.0106. 40% 0.05 0.1703 0.0080 2.2286 0.0574 3.3890 0.0908 2.6948 0.0192. 50% 0.1 0.0014 0.0111 2.0864 0.0519 2.9204 0.0668 2.5322 0.0201. From the above table, we can see that parameters are inaccurate as the white noise level increased. This result means that our iterative estimation method can not be well used. The reason of this consequence is that some of the parameters in 2D-HR type model have the similar solution of x. After the step of noising and smoothing, our solution have a little change which could effect the results of iterative estimation process. In the next figure, we shows that when white noise is 10%, the solution of true value and the solution of estimated value.. 23.

(27) (a) Solution of the true parameter of x. (b) The value of x: plot by the parameter of estimation result. Figure 5.11: On 10% white noise level, the solution of the true parameter and the estimated parameter of x at different time t. 24.

(28) 6. Conclusion. After our work, we present a process to estimate the parameters. Although this paper only discuss the 2D-HR type model, we can expand the iterative process to a general partially observed model. With the consequents in above sections, we have some obstacle in this process. When the observation is without noise, the parameters in partially observed 2DHR model can be well estimated in most of the parameter settings with our iterative estimation process. But when we put noise into the observation, the estimation result will be vary poor. This iterative estimation process may not be vary robust. To this problem, we have 2 ways which may be able to fix it. The first way is reducing the size of ∆t. With more number of data, we can get a smoothing result more similar to the original solution. Use this result, we can estimate the parameters just like we estimated the parameters with an nonnoised observation. The second way is to change the method of smoothing. Not only the kernel regression estimator, but also spline estimator can reduce the noise of our observations. Spline is the main smoothing method in Brunel’s two-way estimator[1]. With different smoothing method, our iterative estimation process might be improved. After all, we wish this process can be used in the real world. With the observations of membrane potential, we can estimate the parameters in 2D-HR type model, and decide which type the neuron is. By this work, we can distinguish the type of a single neuron not only by eyes. Computers can also do this work by parameter estimating. This is our future vision.. 25.

(29) 7 7.1. Some Tables The simulation result while ∆t decreases with k=40, the initial c0 = 10. Table 7.7: Simulation results while ∆t decreases (in the blue area) : 2D-HR model with a = 0, b = 2, c = 3, d = 2.5, T = 50, =0, k = 40, and initial c0 = 10. ∆t. aˆk. bˆk. cˆk. dˆk. bias of cˆk and c. 2−1 2−2 2−3 2−4 2−5 2−6 2−7 2−8 2−9 2−10. -0.1298 -0.0520 -0.0221 -0.0100 -0.0031 -0.0008 -0.0002 -0.0000 -0.0000 -0.0000. 2.8035 1.9489 1.9861 2.0051 2.0023 2.0006 2.0002 2.0000 2.0000 2.0000. 5.3966 2.9855 3.0049 3.0223 3.0085 3.0023 3.0006 3.0001 3.0000 3.0000. 2.6091 2.4148 2.4685 2.4905 2.4974 2.4993 2.4998 2.5000 2.5000 2.5000. 2.3966 -0.0145 0.0049 0.0223 0.0085 0.0023 0.0006 0.0001 0.0000 0.0000. 26.

(30) Table 7.8: Simulation results while ∆t decreases (in the green area) : 2D-HR model with a = −1.1, b = 2, c = 3, d = 2.4, T = 50, =0, k = 40, and initial c0 = 10. ∆t −1. 2 2−2 2−3 2−4 2−5 2−6 2−7 2−8 2−9 2−10 2−11 2−12 2−13 2−14 2−15. aˆk. bˆk. cˆk. dˆk. bias of cˆk and c. -1.1495 -1.1437 -1.1523 -1.1383 -1.1195 -1.1127 -1.1109 -1.1105 -1.1104 -1.1104 -1.1103 -1.1104 -1.1104 -1.1104 -1.1104. 4.8755 4.2499 3.6558 2.8999 2.4076 2.2531 2.2130 2.2031 2.2008 2.2002 2.2001 2.2001 2.2001 2.2001 2.2001. 7.4566 5.8251 4.8692 3.9781 3.4338 3.2658 3.2223 3.2116 3.2091 3.2085 3.2084 3.2084 3.2084 3.2084 3.2084. 3.2878 3.0380 2.8608 2.6436 2.5064 2.4638 2.4528 2.4500 2.4494 2.4493 2.4492 2.4492 2.4492 2.4492 2.4492. 4.4566 2.8251 1.8692 0.9781 0.4338 0.2658 0.2223 0.2116 0.2091 0.2085 0.2084 0.2084 0.2084 0.2084 0.2084. Table 7.9: Simulation results while ∆t decreases (in the red area) : 2D-HR model with a = −1.3, b = 2, c = 3, d = 2.6, T = 50, =0, k = 40, and initial c0 = 10. ∆t. aˆk. bˆk. cˆk. dˆk. bias of cˆk and c. 2−1 2−2 2−3 2−4 2−5 2−6 2−7 2−8 2−9 2−10 2−11 2−12 2−13 2−14 2−15. -1.5097 -1.4577 -1.5643 -1.6838 -1.7313 -1.7436 -1.7471 -1.7483 -1.7488 -1.7490 -1.7491 -1.7492 -1.7492 -1.7492 -1.7492. 6.2225 6.1577 6.9150 7.2555 7.3164 7.3266 7.3310 7.3337 7.3353 7.3361 7.3366 7.3368 7.3369 7.3370 7.3370. 8.7395 8.3006 10.3245 12.3150 12.8089 12.8878 12.9158 12.9320 12.9413 12.9462 12.9488 12.9502 12.9508 12.9512 12.9513. 4.0871 3.9221 4.1158 4.1500 4.1379 4.1346 4.1351 4.1361 4.1368 4.1372 4.1374 4.1375 4.1375 4.1376 4.1376. 5.7395 5.3006 7.3245 9.3150 9.8089 9.8878 9.9158 9.9320 9.9413 9.9462 9.9488 9.9502 9.9508 9.9512 9.9513. 27.

(31) 7.2. The simulation result while ∆t decreases with k=80, the initial c0 = 10. Table 7.10: Simulation results while ∆t decreases (in the blue area) : 2D-HR model with a = 0, b = 2, c = 3, d = 2.5, T = 50, =0, k = 80, and initial c0 = 10. ∆t. aˆk. bˆk. cˆk. dˆk. bias of cˆk and c. 2−1 2−2 2−3 2−4 2−5 2−6 2−7 2−8 2−9 2−10. -0.1298 -0.0520 -0.0221 -0.0100 -0.0031 -0.0008 -0.0002 -0.0000 -0.0000 -0.0000. 2.8035 1.9489 1.9861 2.0051 2.0023 2.0006 2.0002 2.0000 2.0000 2.0000. 5.3966 2.9855 3.0049 3.0223 3.0085 3.0023 3.0006 3.0001 3.0000 3.0000. 2.6091 2.4148 2.4685 2.4905 2.4974 2.4993 2.4998 2.5000 2.5000 2.5000. 2.3966 -0.0145 0.0049 0.0223 0.0085 0.0023 0.0006 0.0001 0.0000 0.0000. 28.

(32) Table 7.11: Simulation results while ∆t decreases (in the green area) : 2D-HR model with a = −1.1, b = 2, c = 3, d = 2.4, T = 50, =0, k = 80, and initial c0 = 10. ∆t −1. 2 2−2 2−3 2−4 2−5 2−6 2−7 2−8 2−9 2−10 2−11 2−12 2−13 2−14 2−15. aˆk. bˆk. cˆk. dˆk. bias of cˆk and c. -1.1495 -1.1437 -1.1522 -1.1366 -1.1132 -1.1044 -1.1020 -1.1014 -1.1012 -1.1012 -1.1012 -1.1012 -1.1012 -1.1012 -1.1012. 4.8755 4.2499 3.6555 2.8670 2.2859 2.0916 2.0405 2.0276 2.0244 2.0236 2.0234 2.0234 2.0233 2.0233 2.0233. 7.4566 5.8251 4.8689 3.9419 3.3061 3.0975 3.0426 3.0288 3.0254 3.0245 3.0243 3.0243 3.0243 3.0243 3.0243. 3.2878 3.0380 2.8607 2.6359 2.4768 2.4242 2.4104 2.4069 2.4060 2.4058 2.4058 2.4058 2.4058 2.4058 2.4058. 4.4566 2.8251 1.8689 0.9419 0.3061 0.0975 0.0426 0.0288 0.0254 0.0245 0.0243 0.0243 0.0243 0.0243 0.0243. Table 7.12: Simulation results while ∆t decreases (in the red area) : 2D-HR model with a = −1.3, b = 2, c = 3, d = 2.6, T = 50, =0, k = 80, and initial c0 = 10. ∆t. aˆk. bˆk. cˆk. dˆk. bias of cˆk and c. 2−1 2−2 2−3 2−4 2−5 2−6 2−7 2−8 2−9 2−10 2−11 2−12 2−13 2−14 2−15. -1.5097 -1.4577 -1.5643 -1.6838 -1.7313 -1.7436 -1.7471 -1.7483 -1.7488 -1.7490 -1.7491 -1.7492 -1.7492 -1.7492 -1.7492. 6.2225 6.1577 6.9150 7.2555 7.3164 7.3266 7.3310 7.3337 7.3353 7.3361 7.3366 7.3368 7.3369 7.3370 7.3370. 8.7395 8.3006 10.3245 12.3150 12.8089 12.8878 12.9158 12.9320 12.9413 12.9462 12.9488 12.9502 12.9508 12.9512 12.9513. 4.0871 3.9221 4.1158 4.1500 4.1379 4.1346 4.1351 4.1361 4.1368 4.1372 4.1374 4.1375 4.1375 4.1376 4.1376. 5.7395 5.3006 7.3245 9.3150 9.8089 9.8878 9.9158 9.9320 9.9413 9.9462 9.9488 9.9502 9.9508 9.9512 9.9513. 29.

(33) 7.3. The simulation result while ∆t decreases with k=40, the initial c0 = 3. Table 7.13: Simulation results while ∆t decreases (in the blue area) : 2D-HR model with a = 0, b = 2, c = 3, d = 2.5, T = 50, =0, k = 40, and initial c0 = 3. ∆t. aˆk. bˆk. cˆk. dˆk. bias of cˆk and c. 2−1 2−2 2−3 2−4 2−5 2−6 2−7 2−8 2−9 2−10. -0.1298 -0.0520 -0.0221 -0.0100 -0.0031 -0.0008 -0.0002 -0.0000 -0.0000 -0.0000. 2.8035 1.9489 1.9861 2.0051 2.0023 2.0006 2.0002 2.0000 2.0000 2.0000. 5.3966 2.9855 3.0049 3.0223 3.0085 3.0023 3.0006 3.0001 3.0000 3.0000. 2.6091 2.4148 2.4685 2.4905 2.4974 2.4993 2.4998 2.5000 2.5000 2.5000. 2.3966 -0.0145 0.0049 0.0223 0.0085 0.0023 0.0006 0.0001 0.0000 0.0000. 30.

(34) Table 7.14: Simulation results while ∆t decreases (in the green area) : 2D-HR model with a = −1.1, b = 2, c = 3, d = 2.4, T = 50, =0, k = 40, and initial c0 = 3. ∆t −1. 2 2−2 2−3 2−4 2−5 2−6 2−7 2−8 2−9 2−10 2−11 2−12 2−13 2−14 2−15. aˆk. bˆk. cˆk. dˆk. bias of cˆk and c. -1.1495 -1.1437 -1.1522 -1.1353 -1.1113 -1.1030 -1.1007 -1.1002 -1.1000 -1.1000 -1.1000 -1.1000 -1.1000 -1.1000 -1.1000. 4.8755 4.2499 3.6552 2.8414 2.2496 2.0637 2.0160 2.0040 2.0010 2.0002 2.0001 2.0000 2.0000 2.0000 2.0000. 7.4566 5.8251 4.8685 3.9139 3.2682 3.0685 3.0172 3.0043 3.0011 3.0003 3.0001 3.0000 3.0000 3.0000 3.0000. 3.2878 3.0380 2.8607 2.6298 2.4679 2.4173 2.4043 2.4011 2.4003 2.4001 2.4000 2.4000 2.4000 2.4000 2.4000. 4.4566 2.8251 1.8685 0.9139 0.2682 0.0685 0.0172 0.0043 0.0011 0.0003 0.0001 0.0000 -0.0000 -0.0000 -0.0000. Table 7.15: Simulation results while ∆t decreases (in the red area) : 2D-HR model with a = −1.3, b = 2, c = 3, d = 2.6, T = 50, =0, k = 40, and initial c0 = 3. ∆t. aˆk. bˆk. cˆk. dˆk. bias of cˆk and c. 2−1 2−2 2−3 2−4 2−5 2−6 2−7 2−8 2−9 2−10 2−11 2−12 2−13 2−14 2−15. -1.5097 -1.4577 -1.5643 -1.6838 -1.7297 -1.5932 -1.3441 -1.3090 -1.3021 -1.3004 -1.3000 -1.2999 -1.2999 -1.2999 -1.2999. 6.2225 6.1577 6.9150 7.2555 7.3176 6.5717 2.8693 2.1813 2.0421 2.0090 2.0007 1.9988 1.9981 1.9982 1.9979. 8.7395 8.3006 10.3245 12.3148 12.7607 9.0707 3.9011 3.1875 3.0436 3.0093 3.0007 2.9987 2.9980 2.9981 2.9979. 4.0871 3.9221 4.1158 4.1500 4.1424 4.1173 2.9136 2.6657 2.6153 2.6032 2.6002 2.5996 2.5993 2.5993 2.5993. 5.7395 5.3006 7.3245 9.3148 9.7607 6.0707 0.9011 0.1875 0.0436 0.0093 0.0007 -0.0013 -0.0020 -0.0019 -0.0021. 31.

(35) 8. Reference. 1 Brunel (2008). Parameter estimation of ODE’s via nonparametric estimators 2 Gugushvili. consistent parameter estimation for systems of ODEs: Bypassing numerical integration via smoothing 3 Nadaraya E. A. (1964). ”On Estimating Regression”. Theory of Probability and its Applications 9 (1): 141V2. 4 Ramsay (2007). Parameter estimation for differential equations a generalized smoothing approach. 32.

(36)

參考文獻

相關文件

The Model-Driven Simulation (MDS) derives performance information based on the application model by analyzing the data flow, working set, cache utilization, work- load, degree

To complete the “plumbing” of associating our vertex data with variables in our shader programs, you need to tell WebGL where in our buffer object to find the vertex data, and

Finally, we use the jump parameters calibrated to the iTraxx market quotes on April 2, 2008 to compare the results of model spreads generated by the analytical method with

Interestingly, the periodicity in the intercept and alpha parameter of our two-stage or five-stage PGARCH(1,1) DGPs does not seem to have any special impacts on the model

Microphone and 600 ohm line conduits shall be mechanically and electrically connected to receptacle boxes and electrically grounded to the audio system ground point.. Lines in

We showed that the BCDM is a unifying model in that conceptual instances could be mapped into instances of five existing bitemporal representational data models: a first normal

(2) We emphasized that our method uses compressed video data to train and detect human behavior, while the proposed method of [19] Alireza Fathi and Greg Mori can only

Besides, we also classify the existing RFID protection mechanisms to solve the different personal privacy threats in our security threat model.. The flowchart of security threat