Question 7.29 (Casella and Berger, 2002)
An alternative to the model in Example 7.2.17 (Casella and Berger, 2002) is the following: we observe \((Y_{i}, X_{i})\) , \(i = 1, 2, \cdots, n\) , where \(Y_{i} \sim \mathcal{P}\mathrm{oisson}(m\beta\tau_{i})\) and \((X_{1}, \cdots, X_{n}) \sim \mathcal{M}\mathrm{ultinomial}(m; \boldsymbol{\tau})\) , with \(\boldsymbol{\tau} = (\tau_{1}, \cdots, \tau_{n})\) such that \(\sum_{i=1}^{n} \tau_{i} = 1\) . Thus, here, for instance, we assume that population counts are multinomial allocations instead of Poisson counts. (Consider \(m = \sum_{i=1}^{n} x_{i}\) as known.)
Show that the joint density of \(\textbf{Y} = (Y_{1}, \cdots, Y_{n})\) and \(\textbf{X} = (X_{1}, \cdots, X_{n})\) is
\[\begin{equation*}
\begin{array}{rclll}
f(\textbf{y}, \textbf{x}\mid \beta, \boldsymbol{\tau}) & = & \displaystyle\prod_{i=1}^{n}\dfrac{e^{-m\beta\tau_{i}}(m\beta\tau_{i})^{y_{i}}}{y_{i}!}m!\dfrac{\tau_{i}^{x_{i}}}{x_{i}!}.
\end{array}
\end{equation*}\]
Solution:
Assuming that \(X_{i}\) and \(Y_{i}\) are mutually independent for all \(i = 1, \cdots, n\) , the joint probability function is given by the product of the probability functions, that is,
\[\begin{equation*}
\begin{array}{rclll}
f(\textbf{y}, \textbf{x}\mid \beta, \boldsymbol{\tau}) & = & \left[ \displaystyle\prod_{i=1}^{n}f(y_{i} \mid \beta, \boldsymbol{\tau})\right] f(\textbf{x}\mid \boldsymbol{\tau})\\[14pt]
& = & \left[\displaystyle\prod_{i=1}^{n}\dfrac{e^{-m\beta\tau_{i}}(m\beta\tau_{i})^{y_{i}}}{y_{i}!}\right]m!\dfrac{\tau_{1}^{x_{1}} \cdots \tau_{n}^{x_{n}}}{x_{1}!\cdots x_{n}!}. \\[14pt]
& = & \displaystyle\prod_{i=1}^{n}\dfrac{e^{-m\beta\tau_{i}}(m\beta\tau_{i})^{y_{i}}}{y_{i}!}m!\dfrac{\tau_{i}^{x_{i}}}{x_{i}!}. \\[8pt]
\end{array}
\end{equation*}\]
If the complete data are observed, show that the maximum likelihood estimators (MLEs) are given by
\[\begin{equation*}
\begin{array}{rclll}
\hat{\beta} & = & \dfrac{\sum_{i=1}^{n}y_{i}}{\sum_{i=1}^{n}x_{i}} \qquad \mbox{e} \qquad \hat{\tau}_{j} & = & \dfrac{x_{j} + y_{j}}{\sum_{i=1}^{n}x_{i} + y_{i}}, \qquad j=1, 2, \cdots, n.
\end{array}
\end{equation*}\]
Solution:
Initially, the log-likelihood function is given by:
\[\begin{equation}
\begin{array}{rclll}
\log [f(\textbf{y}, \textbf{x}\mid \beta, \boldsymbol{\tau})] & = & \log\left[\displaystyle\prod_{i=1}^{n}\dfrac{e^{-m\beta\tau_{i}}(m\beta\tau_{i})^{y_{i}}}{y_{i}!}m!\dfrac{\tau_{i}^{x_{i}}}{x_{i}!}\right] \\[10pt]
& = & \displaystyle\sum_{i=1}^{n}\log\left[\dfrac{e^{-m\beta\tau_{i}}(m\beta\tau_{i})^{y_{i}}}{y_{i}!}m!\dfrac{\tau_{i}^{x_{i}}}{x_{i}!}\right] \\[8pt]
& = & \displaystyle\sum_{i=1}^{n}\left[-m\beta\tau_{i} + \log(m\beta\tau_{i})^{y_{i}}-\log{(y_{i}!)} + \log (m!) + \log(\tau_{i}^{x_{i}}) - \log{(x_{i}!)}\right] \\[8pt]
& = & \displaystyle\sum_{i=1}^{n}\left[-m\beta\tau_{i} + {y_{i}}\log(m\beta\tau_{i})-\log{(y_{i}!)} + \log (m!) + {x_{i}}\log(\tau_{i}) - \log{(x_{i}!)}\right] \\[10pt]
\end{array}
\label{equation1}
\end{equation}\]
Where the first partial derivatives with respect to \(\beta\) and \(\tau_{j}\) are given by the following expressions.
\[\begin{equation*}
\begin{array}{rclll}
\dfrac{\partial}{\partial\beta}\log[f(\textbf{y}, \textbf{x}\mid \beta, \boldsymbol{\tau})] & = & \dfrac{\partial}{\partial\beta}\displaystyle\sum_{i=1}^{n}\left[-m\beta\tau_{i} + {y_{i}}\log(m\beta\tau_{i})-\log{(y_{i}!x_{i}!)} + \log (m!) + {x_{i}}\log(\tau_{i})\right] \\[10pt]
& = & \displaystyle\sum_{i=1}^{n}\dfrac{\partial}{\partial\beta}\left[-m\beta\tau_{i} + {y_{i}}\log(m\beta\tau_{i})-\log{(y_{i}!x_{i}!)} + \log (m!) + {x_{i}}\log(\tau_{i})\right] \\[10pt]
& = & \displaystyle\sum_{i=1}^{n}\left[-m\tau_{i} + {y_{i}}\dfrac{m\tau_{i}}{m\beta\tau_{i}}\right] \\[10pt]
& = & -m\displaystyle\sum_{i=1}^{n}\tau_{i} + \dfrac{1}{\beta}\displaystyle\sum_{i=1}^{n}{y_{i}}\\[12pt]
\end{array}
\end{equation*}\]
Where, as per the initial hypothesis, \(m = \sum_{i=1}^{n} x_{i}\) and \(\sum_{i=1}^{n} \tau_{i} = 1\) , then:
\[\begin{equation*}
\begin{array}{rclll}
\dfrac{\partial}{\partial\beta}\log[f(\textbf{y}, \textbf{x}\mid \beta, \boldsymbol{\tau})]
& = & -\displaystyle\sum_{i=1}^{n}x_{i} + \dfrac{1}{\beta}\displaystyle\sum_{i=1}^{n}{y_{i}}\\[10pt]
\end{array}
\end{equation*}\]
\[\begin{equation*}
\begin{array}{rclll}
\dfrac{\partial}{\partial\tau_{j}}\log[f(y_{j}, x_{j}\mid \beta, \tau_{j})] & = & \dfrac{\partial}{\partial\tau_{j}}\left[-m\beta\tau_{j} + {y_{j}}\log(m\beta\tau_{j})-\log{(y_{j}!x_{j}!)} + \log (m!) + {x_{j}}\log(\tau_{j})\right] \\[6pt]
& = & -m\beta + {y_{j}}\dfrac{m\beta}{m\beta\tau_{j}} + {x_{j}}\dfrac{1}{\tau_{j}} \\[8pt]
& = & -m\beta + {y_{j}}\dfrac{1}{\tau_{j}} + {x_{j}}\dfrac{1}{\tau_{j}} \\[8pt]
& = & -m\beta + \dfrac{y_{j} + x_{j}}{\tau_{j}} \\[8pt]
\end{array}
\end{equation*}\]
Now, by setting these partial derivatives equal to 0, we will obtain as a solution the following maximum likelihood estimators:
\[\begin{equation*}
\begin{array}{rclllll}
\dfrac{\partial}{\partial\beta}\log[f(\textbf{y}, \textbf{x}\mid \beta, \boldsymbol{\tau})] & = & 0 \\[8pt]
-\displaystyle\sum_{i=1}^{n}x_{i} + \dfrac{1}{\hat{\beta}}\displaystyle\sum_{i=1}^{n}{y_{i}} & = & 0\\
\dfrac{1}{\hat{\beta}}\displaystyle\sum_{i=1}^{n}{y_{i}} & = & \displaystyle\sum_{i=1}^{n}x_{i} & \Longrightarrow &
\hat{\beta} & = & \dfrac{\displaystyle\sum_{i=1}^{n}{y_{i}}}{\displaystyle\sum_{i=1}^{n}x_{i}}\\[6pt]
\end{array}
\end{equation*}\]
\[\begin{equation}
\begin{array}{rclll}
\dfrac{\partial}{\partial\tau_{j}}\log[f(y_{j}, x_{j}\mid \beta, \boldsymbol{\tau})] & = & 0 \\[6pt]
-m\hat{\beta} + \dfrac{y_{j} + x_{j}}{\hat{\tau}_{j}} & = & 0\\[8pt]
\dfrac{y_{j} + x_{j}}{\hat{\tau}_{j}} & = & m\hat{\beta}\\[6pt]
\hat{\tau}_{j} & = & \dfrac{y_{j} + x_{j}}{m\hat{\beta}}\\[8pt]
\end{array}
\label{tauemv}
\end{equation}\]
Where, considering this previous result and starting from the initial hypothesis that \(\sum_{i=1}^{n} \tau_{i} = 1\) , we have:
\[\begin{equation*}
\begin{array}{rclll}
%0 & = & -m\beta + \dfrac{y_{j} + x_{j}}{\tau_{j}} \\[8pt]
%\tau_{j} & = & \dfrac{y_{j} + x_{j}}{m\beta} \\[8pt]
\displaystyle\sum_{i=1}^{n}\tau_{i} & = & \displaystyle\sum_{i=1}^{n}\left(\dfrac{y_{i} + x_{i}}{m\beta}\right)\\[8pt]
1 & = & \dfrac{1}{m\beta}\displaystyle\sum_{i=1}^{n}\left(y_{i} + x_{i}\right)\\[8pt]
m\beta & = & \displaystyle\sum_{i=1}^{n}\left(y_{i} + x_{i}\right)\\[8pt]
\end{array}
\end{equation*}\]
Thus, the estimator for \(\hat{\tau}_{j}\) will be given by:
\[\begin{equation*}
\begin{array}{rclll}
\hat{\tau}_{j} & = & \dfrac{y_{j} + x_{j}}{\displaystyle\sum_{i=1}^{n}\left(y_{i} + x_{i}\right)}\\[8pt]
\end{array}
\end{equation*}\]
Q.E.D.
Suppose that \(x_{i}\) is missing. Use the fact that \(X_{1} \sim \mathcal{B}\mathrm{inomial}(m, \tau_{1})\) to calculate the expected log-likelihood of the complete data. Show that the EM sequence is given by
\[\begin{equation*}
\begin{array}{rclll}
\hat{\beta}^{(r+1)} & = & \dfrac{\sum_{i=1}^{n}y_{i}}{m\hat{\tau}_{1}^{(r)} + \sum_{i=2}^{n}x_{i}} \qquad \mbox{e} & \\[18pt]
\qquad \hat{\tau}_{j}^{(r+1)} & = & \dfrac{x_{j} + y_{j}}{m\hat{\tau}_{1}^{(r)} + \sum_{i=2}^{n}x_{i} + \sum_{i=1}^{n}y_{i}}, \qquad j=1, 2, \cdots, n.\\[8pt]
\end{array}
\end{equation*}\]
Solution:
Let the vector \((\textbf{x}_{(-1)}, \textbf{y}) = (y_{1}, (x_{2}, y_{2}), \cdots, (x_{n}, y_{n}))\) , which denotes the incomplete data. The expected log-likelihood of the complete data, using the information provided in the problem that \(X_{1} \sim \mathcal{B}\mathrm{inomial}(m, \tau_{1})\) , will be:
\[\begin{equation*}
\begin{array}{lclll}
\mathbb{E}\left[ \log[f(\textbf{y}, \textbf{x}\mid \beta, \boldsymbol{\tau})] \mid \tau_{1}^{(r)}, (\textbf{x}_{(-1)}, \textbf{y}) \right] \; = \\[14pt]
= \; \displaystyle\sum_{x_{1}=0}^{m} \log\left[\displaystyle\prod_{i=1}^{n}\dfrac{e^{-m\beta\tau_{i}}(m\beta\tau_{i})^{y_{i}}}{y_{i}!}m!\dfrac{\tau_{i}^{x_{i}}}{x_{i}!}\right] {m\choose x_{1}} \left(\tau_{1}^{(r)}\right)^{x_{1}}\left(1-\tau_{1}^{(r)}\right)^{m-x_{1}} \\[18pt]
\end{array}
\end{equation*}\]
\[\begin{equation*}
\begin{array}{lllll}
= & \displaystyle\sum_{x_{1}=0}^{m} \displaystyle\sum_{i=1}^{n}\log\left[\dfrac{e^{-m\beta\tau_{i}}(m\beta\tau_{i})^{y_{i}}}{y_{i}!}m!\dfrac{\tau_{i}^{x_{i}}}{x_{i}!}\right] {m\choose x_{1}} \left(\tau_{1}^{(r)}\right)^{x_{1}}\left(1-\tau_{1}^{(r)}\right)^{m-x_{1}} \\[18pt]
= & \displaystyle\sum_{i=1}^{n}\left[-m\beta\tau_{i} + {y_{i}}[\log(m\beta\tau_{i})]-\log{(y_{i}!)} + \log(m!)\right] \displaystyle\sum_{x_{1}=0}^{m}{m\choose x_{1}} \left(\tau_{1}^{(r)}\right)^{x_{1}}\left(1-\tau_{1}^{(r)}\right)^{m-x_{1}} \; + \\[12pt]
& + \; \displaystyle\sum_{i=2}^{n}\left[{x_{i}}\log(\tau_{i}) -\log{(x_{i}!)}\right]\displaystyle\sum_{x_{1}=0}^{m}{m\choose x_{1}} \left(\tau_{1}^{(r)}\right)^{x_{1}}\left(1-\tau_{1}^{(r)}\right)^{m-x_{1}} \; +
\\[12pt]
& + \displaystyle\sum_{x_{1}=0}^{m}\left[{x_{1}}\log(\tau_{1}) -\log{(x_{1}!)} \right]{m\choose x_{1}} \left(\tau_{1}^{(r)}\right)^{x_{1}}\left(1-\tau_{1}^{(r)}\right)^{m-x_{1}} \\[12pt]
= & \displaystyle\sum_{i=1}^{n}\left[-m\beta\tau_{i} + {y_{i}}[\log(m\beta\tau_{i})]-\log{(y_{i}!)} + \log(m!)\right] \; + \\[12pt]
& + \; \displaystyle\sum_{i=2}^{n}\left[{x_{i}}\log(\tau_{i}) -\log{(x_{i}!)}\right] + \displaystyle\sum_{x_{1}=0}^{m}\left[{x_{1}}\log(\tau_{1}) - \log{(x_{1}!)}\right]{m\choose x_{1}} \left(\tau_{1}^{(r)}\right)^{x_{1}}\left(1-\tau_{1}^{(r)}\right)^{m-x_{1}} \\[12pt]
= & \displaystyle\sum_{i=1}^{n}\left[-m\beta\tau_{i} + {y_{i}}[\log(m\beta\tau_{i})]-\log{(y_{i}!)} + \log(m!)\right] \; + \\[12pt]
& + \; \displaystyle\sum_{i=2}^{n}\left[{x_{i}}\log(\tau_{i}) -\log{(x_{i}!)}\right] + \displaystyle\sum_{x_{1}=0}^{m}\left[{x_{1}}\log(\tau_{1})\right]{m\choose x_{1}} \left(\tau_{1}^{(r)}\right)^{x_{1}}\left(1-\tau_{1}^{(r)}\right)^{m-x_{1}} \; - \\[12pt]
& - \;\displaystyle\sum_{x_{1}=0}^{m} \left[\log{(x_{1}!)} \right]{m\choose x_{1}} \left(\tau_{1}^{(r)}\right)^{x_{1}}\left(1-\tau_{1}^{(r)}\right)^{m-x_{1}} \\[12pt]
= & \displaystyle\sum_{i=1}^{n}\left[-m\beta\tau_{i} + {y_{i}}[\log(m) + \log(\beta) + \log(\tau_{i})]-\log{(y_{i}!)} + \log(m!)\right] \; + \\[12pt]
& + \; \displaystyle\sum_{i=2}^{n}\left[{x_{i}}\log(\tau_{i}) -\log{(x_{i}!)}\right] + \log(\tau_{1})\displaystyle\sum_{x_{1}=0}^{m}\left[{x_{1}}\right]{m\choose x_{1}} \left(\tau_{1}^{(r)}\right)^{x_{1}}\left(1-\tau_{1}^{(r)}\right)^{m-x_{1}} \; - \\[12pt]
& - \;\displaystyle\sum_{x_{1}=0}^{m} \left[\log{(x_{1}!)} \right]{m\choose x_{1}} \left(\tau_{1}^{(r)}\right)^{x_{1}}\left(1-\tau_{1}^{(r)}\right)^{m-x_{1}} \\[12pt]
= & \left( \displaystyle\sum_{i=1}^{n}\left[-m\beta\tau_{i} + {y_{i}}[\log(m) + \log(\beta) + \log(\tau_{i})]-\log{(y_{i}!)} + \log(m!)\right] \; + \right. \\[12pt]
& + \;\left. \displaystyle\sum_{i=2}^{n}\left[{x_{i}}\log(\tau_{i}) -\log{(x_{i}!)}\right] + \log(\tau_{1})\mathbb{E}\left( X_{1}\mid \tau_{1}^{(r)}\right)\right) \; - \\[12pt]
& - \;\left( \displaystyle\sum_{x_{1}=0}^{m} \left[\log{(x_{1}!)} \right]{m\choose x_{1}} \left(\tau_{1}^{(r)}\right)^{x_{1}}\left(1-\tau_{1}^{(r)}\right)^{m-x_{1}} \right) \\[12pt]
\end{array}
\end{equation*}\]
Since we are calculating this expected value of the log-likelihood with the purpose of maximizing with respect to \(\beta\) and \(\boldsymbol{\tau} = (\tau_{1}, \cdots, \tau_{n})\) , we can ignore the terms in the second set of parentheses. Therefore, we only need to maximize the terms in the first set of parentheses, where we observe that the expected log-likelihood of the complete data is equivalent to the original log-likelihood of the complete data found in part (b), with the exception that \(x_{1}\) is replaced by \(\mathbb{E}\left( X_{1} \mid \tau_{1}^{(r)} \right) = m\tau_{1}^{(r)}\) , as per the initial hypothesis. Thus,
\[\begin{equation*}
\begin{array}{lclll}
\mathbb{E}\left[ \log[f(\textbf{y}, \textbf{x}\mid \beta, \boldsymbol{\tau})] \mid \tau_{1}^{(r)}, (\textbf{x}_{(-1)}, \textbf{y}) \right] \; = \\[10pt]
= \; \left( \displaystyle\sum_{i=1}^{n}\left[-m\beta\tau_{i} + {y_{i}}[\log(m) + \log(\beta) + \log(\tau_{i})]-\log{(y_{i}!)} + \log(m!)\right] \; + \right. \\[12pt]
\; + \;\left. \displaystyle\sum_{i=2}^{n}\left[{x_{i}}\log(\tau_{i}) -\log{(x_{i}!)}\right] + \log(\tau_{1})m\tau_{1}^{(r)}\right) - \left( \displaystyle\sum_{x_{1}=0}^{m} \left[\log{(x_{1}!)} \right]{m\choose x_{1}} \left(\tau_{1}^{(r)}\right)^{x_{1}}\left(1-\tau_{1}^{(r)}\right)^{m-x_{1}} \right) \\[12pt]
\end{array}
\end{equation*}\]
Therefore, in the \(r\) -th iteration, the maximum likelihood estimators (MLEs) are simply a variation of the MLEs found in part (b), replacing \(x_{1}\) by \(\mathbb{E}\left( X_{1} \mid \tau_{1}^{(r)} \right) = m\tau_{1}^{(r)}\) , that is,
\[\begin{equation*}
\begin{array}{rclllll}
\hat{\beta}^{(r+1)} & = & \dfrac{\sum_{i=1}^{n}y_{i}}{x_{1}+ \sum_{i=2}^{n}x_{i}} & = & \dfrac{\sum_{i=1}^{n}y_{i}}{m\hat{\tau}_{1}^{(r)} + \sum_{i=2}^{n}x_{i}} \\[18pt]
\qquad \hat{\tau}_{j}^{(r+1)} & = & \dfrac{x_{j} + y_{j}}{x_{1} + \sum_{i=2}^{n}x_{i} + \sum_{i=1}^{n}y_{i}} & = & \dfrac{x_{j} + y_{j}}{m\hat{\tau}_{1}^{(r)} + \sum_{i=2}^{n}x_{i} + \sum_{i=1}^{n}y_{i}}\\[10pt]
\end{array}
\end{equation*}\]
Thus, we have that the EM sequence is given by:
\[\begin{equation*}
\begin{array}{rclllll}
\hat{\beta}^{(r+1)} & = & \dfrac{\sum_{i=1}^{n}y_{i}}{m\hat{\tau}_{1}^{(r)} + \sum_{i=2}^{n}x_{i}}, & \qquad \hat{\tau}_{1}^{(r+1)} & = & \dfrac{m\hat{\tau}_{1}^{(r)} + y_{1}}{m\hat{\tau}_{1}^{(r)} + \sum_{i=2}^{n}x_{i} + \sum_{i=1}^{n}y_{i}} \\[18pt]
\end{array}
\end{equation*}\] \[\begin{equation*}
\begin{array}{rclllll}
\qquad \hat{\tau}_{j}^{(r+1)} & = & \dfrac{x_{j} + y_{j}}{m\hat{\tau}_{1}^{(r)} + \sum_{i=2}^{n}x_{i} + \sum_{i=1}^{n}y_{i}}, \qquad j=2, \cdots, n.\\[15pt]
\end{array}
\end{equation*}\]
Q.E.D.
Use this model to find the maximum likelihood estimators (MLEs) for the data from Exercise 7.28 (Casella and Berger, 2002), first assuming that you have all the data and then assuming that \(x_{1} = 3540\) is missing.
Solution:
To find the maximum likelihood estimators (MLEs) for \(\hat{\beta}\) and \(\boldsymbol{\hat{\tau}}\) in this question, we will use the programming language (R Core Team, 2023) to implement the algorithm.
Thus, the MLEs found using the complete data were the following: \[\begin{equation*}
\begin{array}{rclll}
\hat{\beta} & = & 0.0008413892 \\
\hat{\boldsymbol{\tau}} & = & (0.06337310, 0.06374873, 0.06689681, 0.04981487, 0.04604075, 0.04883109, \\
&& 0.07072460, 0.01776164, 0.03416388, 0.01695673, 0.02098127, 0.01878119, \\
&& 0.05621836, 0.09818091, 0.09945087, 0.05267677, 0.08896918, 0.08642925)\\
\end{array}
\end{equation*}\]
Now, the MLEs found using the incomplete data, and considering \(m = \sum_{i} x_{i}\) , were:
\[\begin{equation*}
\begin{array}{rclll}
\hat{\beta} & = & 0.0008415534 \\
\hat{\boldsymbol{\tau}} & = & (0.06319044, 0.06376116, 0.06690986, 0.04982459, 0.04604973, 0.04884062, \\
&& 0.07073839, 0.01776510, 0.03417054, 0.01696004, 0.02098536, 0.01878485, \\
&& 0.05622933, 0.09820005, 0.09947027, 0.05268704, 0.08898653, 0.08644610)\\
\end{array}
\end{equation*}\]
It is worth noting that in the code given below, the value of \(\tau_{1}\) used in the second algorithm for the incomplete data is initialized with the \(\tau_{1}\) found in the result of the first algorithm with the complete data. The code is as follows:
# emvs para os dados completos xdata <- c ( 3540 , 3560 , 3739 , 2784 , 2571 , 2729 , 3952 , 993 , 1908 , 948 , 1172 ,
1047 , 3138 , 5485 , 5554 , 2943 , 4969 , 4828
) ydata <- c ( 3 , 4 , 1 , 1 , 3 , 1 , 2 , 0 , 2 , 0 , 1 , 3 , 5 , 4 , 6 , 2 , 5 , 4 ) for ( j in 1 : 500 ) { beta <- sum ( ydata ) / sum ( xdata )
tau <- c ( ( xdata + ydata ) / ( sum ( xdata ) + sum ( ydata ) ) )
} beta
[1] 0.06337310 0.06374873 0.06689681 0.04981487 0.04604075 0.04883109
[7] 0.07072460 0.01776164 0.03416388 0.01695673 0.02098127 0.01878119
[13] 0.05621836 0.09818091 0.09945087 0.05267677 0.08896918 0.08642925
# emvs para os dados imcompletos xdatam <- c ( 3560 , 3739 , 2784 , 2571 , 2729 , 3952 , 993 , 1908 , 948 , 1172 ,
1047 , 3138 , 5485 , 5554 , 2943 , 4969 , 4828
) ydata <- c ( 3 , 4 , 1 , 1 , 3 , 1 , 2 , 0 , 2 , 0 , 1 , 3 , 5 , 4 , 6 , 2 , 5 , 4 ) xdata <- c ( mean ( xdatam ) , xdatam ) i <- 0 for ( j in 1 : 500 ) { xdata <- c ( sum ( xdata ) * tau [ 1 ] , xdatam )
beta <- sum ( ydata ) / sum ( xdata )
tau <- c ( ( xdata + ydata ) / ( sum ( xdata ) + sum ( ydata ) ) )
i <- i + 1
} beta
[1] 0.06319044 0.06376116 0.06690986 0.04982459 0.04604973 0.04884062
[7] 0.07073839 0.01776510 0.03417054 0.01696004 0.02098536 0.01878485
[13] 0.05622933 0.09820005 0.09947027 0.05268704 0.08898653 0.08644610