কোর্স:পাসপ্র/সমস্যা ৩.৩.১

ডিপথট থেকে
<< সমস্যা ৩.২.১ কোর্স:পাসপ্র সমস্যা ৩.৩.২ >>
প্রশিক্ষক: সমস্যা ৩.৩.১

সমস্যা

অনেক সময়ই আমরা যে মডেলটি ফিট করতে চাচ্ছি তার রাশিগুলো রৈখিক হয় না। যেমন, নিচের সমীকরণের মতো যদি আমরা একটি সাইনুসয়েডের দশা নির্ণয় করতে চাই।

<math>x[n] = \sin(2\pi f n+\phi) + w[n]</math>

যেখানে <math>n=0,1,2,...,N-1</math> এবং <math>w[n] \rightarrow \mathcal{N}(0,\sigma^2)</math>

  1. টেইলর সম্প্রসারণের মাধ্যমে একে <math>\phi_0</math> এর চারপাশে একটি রৈখিক সমীকরণে রূপান্তরিত করো। এই মডেলের <math>\phi</math> রাশির জন্য ন্যূভেনি এবং সিআরএলবি নির্ণয় করো। এই রৈখিক নিরূপকের টাল নির্ণয় করো। কোন শর্তে এটি একটি নিটাল নিরূপক হবে?
  2. এখন উপরের সমস্যাটি কম্পিউটার প্রোগ্রামের মাধ্যমে সমাধান করো। প্রথমে <math>N=500, \sigma=0.5, f=0.2, \phi=1</math> ধরে একটি উপাত্ত তৈরি করো। এরপর <math>\phi_0=[0.5, 0.8, 1.2]</math> ধরে রৈখিক মডেল ফিটিং এর মাধ্যমে <math>\phi</math> নিরূপণ করো। দেখো <math>\phi_0</math> কোন মানের জন্য নিরূপিত মান প্রকৃত মানের সবচেয়ে কাছাকাছি হয়।
  3. এই পরীক্ষাটি ১০০০ বার পুনরাবৃত্তি করে নিরূপিত মানগুলোর হিস্টোগ্রাম তৈরি করো। হিস্টোগ্রামের আকৃতির ব্যাখ্যা করো।

সমাধান ১

প্রথমেই টেইলর সম্প্রসারণের সূত্রটি জেনে নেয়া যাক। যেকোনো অপেক্ষক <math>f(x)</math> কে তার যেকোনো মান <math>a</math> কেন্দ্র করে টেইলর সম্প্রসারণ করলে পাওয়া যাবে:

<math>f(x) = f(a)+\frac {f'(a)}{1!} (x-a)+ \frac{f(a)}{2!} (x-a)^2+\frac{f^{(3)}(a)}{3!}(x-a)^3+ \cdots. = \sum_{n=0} ^ {\infty} \frac {f^{(n)}(a)}{n!} \, (x-a)^{n} </math>

এই সূতে মেনে আমাদের ফাংশনটির টেইলর সম্প্রসারণ করে ১ম ঘাতের পর সবকিছু অগ্রাহ্য করলে পাওয়া যায়:

<math>x[n] \approx \sin(2\pi f n+\phi_0) + (\phi-\phi_0) \cos(2\pi f n+\phi_0) + w[n]</math>

ধরি,

<math>\Delta\phi = \phi-\phi_0</math>
<math>x'[n] = x[n]-\sin(2\pi f n+\phi_0)</math>

তাহলে,

<math>x'[n] = \cos(2\pi f n+\phi_0) \Delta\phi + w[n]</math>

এটি একটি রৈখিক মডেল যা নিচের সমীকরণের মতো:

<math>\mathbf{X}' = \mathbf{H}\boldsymbol{\theta} + \mathbf{W}</math>

যেখানে:

<math>\mathbf{X}' = \begin{bmatrix}

x'[o] \\ x'[1] \\ \vdots \\ x'[N-1] \end{bmatrix} \ \ \mathbf{H} = \begin{bmatrix} \cos(2\pi f.0+\phi_0) \\ \cos(2\pi f.1+\phi_0) \\ \vdots \\ \cos(2\pi f.(N-1)+\phi_0) \end{bmatrix} </math> এবং <math>\boldsymbol{\theta} = \Delta\phi\</math>.

এখন এই রৈখিক মডেলের জন্য একটি ন্যূভেনি নির্ণয় করতে হবে। রৈখিক মডেলের কোনো অজানা রাশির মান নিরূপণের জন্য নিরূপক কেমন হয় সেটা আমরা রৈখিক মডেল অংশে আলোচনা করেছিলাম। সেই অনুযায়ী আমাদের উপরের মডেলের জন্য ন্যূভেনি হবে:

<math>\hat{\boldsymbol{\theta}} = (\mathbf{H}^T\mathbf{H})^{-1} \mathbf{H}^T \mathbf{X}'</math>

যেখান থেকে আমরা লিখতে পারি,

<math>\hat{\phi} = \phi_0 + (\mathbf{H}^T\mathbf{H})^{-1} \mathbf{H}^T \mathbf{X}'</math>

যেহেতু এটি ন্যূভেনি, সেহেতু এর ভেদাঙ্ক সম্ভাব্য সর্বনিম্ন, অর্থাৎ এর ভেদাঙ্কই ক্রামার-রাও নিম্নসীমা। নিরূপক <math>\hat{\boldsymbol{\theta}}</math> যেহেতু একটি মেট্রিক্স সেহেতু এর ভেদাঙ্ক নয়, বরং সহভেদাঙ্ক নির্ণয় করতে হবে। একটি রৈখিক নিরূপকের সহভেদাঙ্কের সমীকরণ হচ্ছে:

<math>C_{\hat{\boldsymbol{\theta}}} = \sigma^2 (\mathbf{H}^T\mathbf{H})^{-1} </math>
<math>\mathbf{H}^T\mathbf{H} = \sum\limits_{n=0}^{N-1} \cos^2(2\pi f n+\phi_0) \approx \frac{N}{2}</math>
<math>\therefore C_{\hat{\boldsymbol{\theta}}} = \frac{2\sigma^2}{N}</math>

তাহলে <math>\phi</math> এর সাপেক্ষে সিআরএলবি হবে <math>2\sigma^2/N</math>

রৈখিক মডেলের জন্য নিরূপক ছিল:

<math>\hat{\boldsymbol{\phi}} = \phi_0 + (\mathbf{H}^T\mathbf{H})^{-1} \mathbf{H}^T \mathbf{X}'</math>

এই নিরূপকের প্রত্যাশা নির্ণয় করে দেখতে হবে সেটি মূল মানের সমান কি না। সমান হলে নিরূপকটি নিটাল, আর সমান না হলে দেখতে হবে কোন শর্তে সেটিকে সমান করা যায়।

<math>\begin{align}

E(\hat{\boldsymbol{\phi}}) &= \phi_0 + (\mathbf{H}^T\mathbf{H})^{-1} E(\mathbf{H}^T \mathbf{X}') \\ &\approx \phi_0 + \frac{2}{N} E\left( \sum\limits_{n=0}^{N-1} \cos(2\pi fn+\phi_0) [\sin(2\pi fn+\phi) + w[n]] \right) \\ &= \phi_0 + \frac{2}{N} \sum\limits_{n=0}^{N-1} \cos(2\pi fn+\phi_0) \sin(2\pi fn+\phi) + \frac{2}{N} \sum\limits_{n=0}^{N-1} E\left(w[n] \cos(2\pi fn+\phi_0) \right) \\ &= \phi_0 + \frac{2}{N} \sum\limits_{n=0}^{N-1} \cos(2\pi fn+\phi_0) \sin(2\pi fn+\phi) \\ &\approx \phi_0 + \frac{1}{N} \sum\limits_{n=0}^{N-1} \sin(4\pi fn+\phi+\phi_0) + \frac{1}{N} \sum\limits_{n=0}^{N-1} \sin(\phi-\phi_0) \\ &= \phi_0 + \frac{1}{N} \sum\limits_{n=0}^{N-1} \sin(\phi-\phi_0) \end{align}</math> এটাই আমাদের নিরূপিত মান। এটি নিশ্চিত প্রকৃত রাশির (<math>\phi</math>) সমান নয়, কারণ এটা আমাদের পূর্বানুমিত মানের (<math>\phi_0</math>) উপর নির্ভর করছে। এখন আমাদের পূর্বানুমিত মান যদি প্রকৃত মানের খুব কাছাকাছি হয় তাহলে <math>\phi-\phi_0</math> খুব ছোটো হবে, এবং আমাদের নিরূপকটি নিটাল হয়ে যাবে:

<math>E(\hat{\boldsymbol{\phi}}) \approx \phi_0 + \frac{1}{N} N (\phi-\phi_0) = \phi_0 + (\phi-\phi_0) \approx \phi</math>

সমাধান ২

নিচের ম্যাটল্যাব প্রোগ্রামটি দিয়ে দশার মান নিরূপণ করা সম্ভব:

clear all; close all; clc;
% data
N=500; sig=0.5; f=0.2; phi=1;
n=[0:(N-1)]';
x=sin(2*pi*f*n+phi);
w=randn(size(x))*sig;
x=x+w;

% linear estimator for phi
phi0=[0.5 0.8 1.2];
phi_est = zeros(size(phi0));
for i=1:length(phi0)
    H=cos(2*pi*f*n+phi0(i));
    xp=x-sin(2*pi*f*n+phi0(i));
    phi_est(i)=phi0(i)+inv(H'*H)*(H'*xp);
end
disp(['Initial values used were (radians) ' num2str(phi0)]);
disp(['Estimated values are (radians) ' num2str(phi_est)]);

পূর্বানুমিত দশা ০.৫, ০.৮ এবং ১.২ রেডিয়ান হলে নিরূপত দশা হয় যথাক্রমে ০.৯৬৭৮, ০.৯৯১৭, এবং ১.০০১৪। দেখা যাচ্ছে পূর্বানুমিত মান প্রকৃত মানের (১) যত কাছাকাছি হয় নিরূপিত মানও ততই প্রকৃতটির কাছাকাছি হয়।

সমাধান ৩

দৈব নয়েজবিশিষ্ট উপাত্তের একাধিক (আসলে অসংখ্য) নমুনা তৈরি করে প্রতিটি নমুনার জন্য একই পরীক্ষা পুনরাবৃত্তি করার প্রক্রিয়াকে মন্টি-কার্লো পদ্ধতি বলা হয়। ঠিক এটাই আমরা এখানে করব। এক্ষেত্রে আমাদের উপের প্রোগ্রামটিতে শুধু একটি পরিবর্তন আসবে: আমরা উপের প্রোগ্রামটি ১০০০ বার রান করাব এবং প্রতি বার গাউসীয় দৈব নয়েজের একটি ভিন্ন স্বাধীন নমুনা নেব। তারপর হিস্টোগ্রাম তৈরি করব। নিচের ম্যাটল্যাব কোড দিয়ে পুরোটা করা সম্ভব:

clear all; close all; clc;
% data
N=500; sig=0.5; f=0.2; phi=1;
n=[0:(N-1)]';
x=sin(2*pi*f*n+phi);
Nrealz = 1000; % no. of Monte-Carlo realizations

% linear estimator for phi
phi0=[0.5 0.8 1.2];
phi_est = zeros(Nrealz,length(phi0));
for i=1:length(phi0)
    H=cos(2*pi*f*n+phi0(i));
    for j=1:Nrealz
        xp=x+sig*randn(size(x))-sin(2*pi*f*n+phi0(i));
        phi_est(j,i)=phi0(i)+inv(H'*H)*(H'*xp);
    end
end
figure;
subplot(131); hist(phi_est(:,1),40);
xlabel('Estimated value of phase (radians)');
ylabel('No. of occurences');
subplot(132); hist(phi_est(:,2),40);
xlabel('Estimated value of phase (radians)');
ylabel('No. of occurences');
subplot(133); hist(phi_est(:,3),40);
xlabel('Estimated value of phase (radians)');
ylabel('No. of occurences');
axes('position',[0,0,1,1],'visible','off');
text(.82,.90, ['Mean = ' num2str(mean(phi_est(3)))] );
text(.82,.87, ['Variance = ' num2str(mean(phi_est(3)))] );
text(.82,.83, '\phi_0 = 1.2' );
text(.54,.90, ['Mean = ' num2str(mean(phi_est(2)))] );
text(.54,.87, ['Variance = ' num2str(mean(phi_est(2)))] );
text(.54,.83, '\phi_0 = 0.8' );
text(.26,.90, ['Mean = ' num2str(mean(phi_est(1)))] );
text(.26,.87, ['Variance = ' num2str(mean(phi_est(1)))] );
text(.26,.83, '\phi_0 = 0.5' );

উপরের স্ক্রিপ্টটি রান করালে নিচের ছবিটি পাওয়া যাবে:

SSP hw423.png

এক্ষেত্রেও দেখা যাচ্ছে পূর্বানুমিত মান প্রকৃত মানের যত কাছাকাছি হয় মন্টি-কার্লো রিয়েলাইজেশন সমূহের গড়ও তত প্রকৃত মানের কাছাকাছি হয়, এবং ভেদাঙ্ক তত কম হয়। উল্লেখ্য ভেদাঙ্ক <math>10^{-3} \ \text{rad}^2</math> এর এককে দেখানো হয়েছে। গড়ের একক রেডিয়ান।

সমস্যা ১ এর সূত্র (<math>\text{CRLB}=2\sigma^2/N</math>) থেকে হিসাব করলে দেখা যাবে <math>\sigma=0.5</math> এর জন্য সিআরএলবি <math>10^{-3}</math> যা আমাদের নিরূপিত ভেদাঙ্কের মানের খুব কাছাকাছি।