使用蒙特卡羅模擬計(jì)算功率，*1部分: 基礎(chǔ)知識(shí)

時(shí)間：2019-12-20作者：北京天演融智軟件有限公司瀏覽：180

功率和樣本大小的計(jì)算是科學(xué)研究計(jì)劃的重要組成部分。可以使用Stata的power命令來計(jì)算許多常用統(tǒng)計(jì)測試的功率和樣本大小需求。但對于較復(fù)雜的模型是沒有簡單公式的，如多層/縱向模型和結(jié)構(gòu)方程模型(SEMs)。蒙特卡羅模擬是計(jì)算復(fù)雜模型的功率和樣本大小要求的一種方法，Stata提供了執(zhí)行此操作所需要的所有工具。甚至可以將模擬集成到Stata的power命令中，可以輕松地為一系列參數(shù)值創(chuàng)建自定義表格和圖形。




比如，下面的自定義程序power simmixed模擬縱向模型的功率，假設(shè)參與者的數(shù)量（級別2）為100到500，每次增加100，每個(gè)參與者有5到6個(gè)觀察值(級別1)。power simmixed還創(chuàng)建了一個(gè)表(未顯示)和下圖來顯示模擬的結(jié)果。




我的同事和我寫了一系列文章展示如何做到這一點(diǎn)。在今天的文章中，將介紹使用模擬計(jì)算功率和樣本大小要求所需的基本工具。在*二篇文章中，將展示如何將模擬集成到Stata的power命令中。 然后，我們將展示線性回歸，邏輯回歸，多級/縱向模型和結(jié)構(gòu)方程模型的具體示例。

基本思路

統(tǒng)計(jì)功率是當(dāng)零假設(shè)為假時(shí)拒絕零假設(shè)的概率。功率的計(jì)算是基于一組假設(shè)，例如樣本大小，alpha水平和特定的替代假設(shè)。例如，我們可能希望計(jì)算t檢驗(yàn)的功率，假設(shè)零假設(shè)的樣本均值為70，替代假設(shè)為75，樣本大小為100，alpha水平為0.05。

使用蒙特卡羅模擬計(jì)算功率的基本步驟是

1. 生成假設(shè)替代假設(shè)為真(例如，均值=75)的數(shù)據(jù)集。

2. 使用數(shù)據(jù)集測試零假設(shè)（例如，測試均值= 70）。

3. 保存測試結(jié)果（例如，“拒絕”或“未拒絕”）。

4. 多次重復(fù)步驟1-3（通常為1,000或更多）。

零假設(shè)被拒絕的比例是我們對統(tǒng)計(jì)能力的估算。在上面的例子中，我們可能在1000次迭代中觀察到834次“拒絕”，這使我們的估算功率為0.834或83.4％。

要執(zhí)行這些步驟，需要熟悉Stata的一些編程工具。以下是本文中介紹的主題列表。 如果您熟悉其中一些主題，可以單擊下面的鏈接跳到不熟悉的主題。

主題列表

標(biāo)量和局部宏

創(chuàng)建偽隨機(jī)數(shù)據(jù)集

存儲(chǔ)模型輸出

如何使用程序創(chuàng)建簡單程序

如何使用程序創(chuàng)建有用的程序

如何使用模擬多次運(yùn)行程序

用power onemean來檢查結(jié)果

標(biāo)量和局部宏

標(biāo)量和本地宏是模擬的重要工具，因?yàn)樗鼈冊试S將數(shù)字臨時(shí)存儲(chǔ)在內(nèi)存中。例如，可以通過輸入將數(shù)字1存儲(chǔ)到名為i的標(biāo)量中




稍后通過輸入來引用此標(biāo)量







如果對標(biāo)量和變量使用相同的名稱，Stata將會(huì)變得混亂?？梢酝ㄟ^為標(biāo)量使用一的名稱來避免這種混淆，或者可以使用scalar()函數(shù)來引用標(biāo)量。




還可以使用本地宏存儲(chǔ)數(shù)字。 例如，可以通過輸入將數(shù)字1存儲(chǔ)到名為i的本地宏。







然后，可以通過在其**引用左單引號(hào)（通常位于“1”鍵左側(cè)的鍵盤上）并使用右單引號(hào)（通常位于Enter鍵左側(cè)）來引用本地宏。）。




本地宏通常用于定義模擬的輸入?yún)?shù)，許多Stata命令的結(jié)果存儲(chǔ)為標(biāo)量。

創(chuàng)建偽隨機(jī)數(shù)據(jù)集

您還需要生成隨機(jī)數(shù)來進(jìn)行模擬。在這里，將展示一些常用的隨機(jī)數(shù)函數(shù)，并參考Stata Functions Reference Manual以獲取完整列表。首先清除Stata的內(nèi)存，將隨機(jī)數(shù)種子設(shè)置為15，然后使用set obs告訴Stata我們想要?jiǎng)?chuàng)建一個(gè)包含200個(gè)觀察值的數(shù)據(jù)集。




可以使用runiform()函數(shù)在區(qū)間(0,1)內(nèi)生成均勻分布的隨機(jī)數(shù)。




可以使用類似的runiformint(a,b)函數(shù)在區(qū)間[a,b]上生成均勻分布的隨機(jī)整數(shù)變量。例如，可以輸入以下命令，在區(qū)間[18,65]內(nèi)生成一個(gè)人的年齡的隨機(jī)值:




可以使用rbinomial（n，p）函數(shù)生成二項(xiàng)式(n,p)隨機(jī)變量，其中n是試驗(yàn)次數(shù)，p是成功概率。例如，可以通過輸入以下命令為女性生成隨機(jī)指示器：







還可以使用rnormal(m,s)函數(shù)從正常密度生成隨機(jī)值，其中均值等于m，標(biāo)準(zhǔn)差等于s。 例如，可以使用以下命令生成重量和高度的變量。在這里，*了72千克的平均值和15的標(biāo)準(zhǔn)偏差，以及170厘米的平均值和10的標(biāo)準(zhǔn)偏差。




在這個(gè)例子中，我們獨(dú)立地生成了weight和height。但是像weight和height這樣的變量很可能是相關(guān)的，可以使用drawnorm來生成相關(guān)的變量。

在下面的示例中，平均值存儲(chǔ)在矩陣m中，標(biāo)準(zhǔn)偏差存儲(chǔ)在矩陣s中，矩陣C中存儲(chǔ)變量之間的相關(guān)性。然后，可以將這些矩陣作為參數(shù)包含在drawingorm選項(xiàng)中以創(chuàng)建相關(guān)變量height和weight.。







上述均值和標(biāo)準(zhǔn)差的估算值與我們的輸入?yún)?shù)相似，下面的估算相關(guān)系數(shù)為0.5049，接近我們在上面的相關(guān)矩陣中*的0.5的值。




存儲(chǔ)模型輸出

許多Stata命令將其結(jié)果存儲(chǔ)為標(biāo)量，宏和矩陣。運(yùn)行命令后，可以通過輸入return list查看存儲(chǔ)結(jié)果的列表。還可以在估算命令（如regress）后輸入ereturn list。

在下面的例子中，我使用ttest命令來測試平均權(quán)重等于70的零假設(shè)。







輸入return list顯示Stata存儲(chǔ)在內(nèi)存中的標(biāo)量列表。




上述均值和標(biāo)準(zhǔn)差的估算值與我們的輸入?yún)?shù)相似，下面的估算相關(guān)系數(shù)為0.5049，接近我們在上面的相關(guān)矩陣中*的0.5的值。







存儲(chǔ)輸出模型

許多Stata命令將其結(jié)果存儲(chǔ)為標(biāo)量，宏和矩陣。運(yùn)行命令后，可以通過輸入return list查看存儲(chǔ)結(jié)果的列表。還可以在估算命令（如regress）后鍵入ereturn list。

在下面的示例中，我使用ttest命令來測試平均權(quán)重等于70的零假設(shè)。




輸入return list會(huì)顯示Stata存儲(chǔ)在內(nèi)存中的標(biāo)量列表。







可以將這些標(biāo)量中的任意一個(gè)存儲(chǔ)到另一個(gè)標(biāo)量中。例如，可以通過輸入將存儲(chǔ)在r(p)中的雙面p值存儲(chǔ)到名為pvalue的標(biāo)量中




如何使用程序創(chuàng)建一個(gè)簡單的程序

您將運(yùn)行生成隨機(jī)數(shù)據(jù)的命令，并在運(yùn)行模擬時(shí)多次測試零假設(shè)。還需要為模擬定義輸入?yún)?shù)并返回假設(shè)檢驗(yàn)的結(jié)果。有效地完成這些任務(wù)的一種方法是使用program定義自己的Stata程序。下面的代碼塊定義了一個(gè)名為myprogram的程序，它接受輸入?yún)?shù)n()，顯示n的值，并返回n的值。




讓我們逐行考慮這個(gè)代碼塊。

**行是capture program drop myprogram。定義程序后，必須先使用program drop將程序從內(nèi)存中刪除，然后才能修改和重新定義程序。因?yàn)楹芸赡茉谕瓿芍皶?huì)多次修改此程序，所以在重新定義程序之前，它將節(jié)省運(yùn)行program drop的時(shí)間。**次運(yùn)行此代碼塊時(shí)不會(huì)定義程序，所以program drop將返回一個(gè)錯(cuò)誤。在program drop之前輸入capture將捕獲此錯(cuò)誤并允許代碼塊繼續(xù)運(yùn)行。

*二行是程序program myprogram, rclass。 開始定義程序myprogram。 定義程序后，可以輸入myprogram，Stata將在program和end之間運(yùn)行所有命令。選項(xiàng)rclass告訴Stata我們想要使用return返回程序中的值。

*三行，15.1版本告訴Stata希望使用Stata 15.1中編寫的功能來運(yùn)行程序。 可以在Stata Programming Reference Manual中了解有關(guān)Stata版本控制的更多信息。

*四行是syntax, n(integer)。 定義了程序的語法。該程序要求用戶輸入逗號(hào)，后跟n()的整數(shù)。n()的值存儲(chǔ)在名為n的本地宏中，可以將此本地宏稱為“n”。

*五行顯示本地宏n的值。

*六行是return scalar N = `n’。此行指示程序以名為N的標(biāo)量形式返回n的值。
最后一行是end。這告訴Stata你已經(jīng)完成了程序myprogram的定義。

讓我們運(yùn)行myprogram，看看它做了什么。




為n()輸入值50，Stata顯示結(jié)果n = 50。
可以輸入return list，并看到myprogram以標(biāo)量r（N）的形式返回n的值。




非常有用！

如何使用程序創(chuàng)建有用的程序

讓我們定義一個(gè)名為simttest的程序，它根據(jù)輸入?yún)?shù)生成一個(gè)隨機(jī)數(shù)據(jù)集，測試零假設(shè)，并返回我們的假設(shè)檢驗(yàn)結(jié)果。 下面的代碼塊僅使用注釋來定義程序。




現(xiàn)在，讓我們在下面的代碼塊中填寫詳細(xì)信息。




syntax的定義包括一個(gè)必需參數(shù)，n()和四個(gè)可選參數(shù)，alpha(),m0(),ma(),和sd()。 可選參數(shù)包含在方括號(hào)中。所有輸入?yún)?shù)都包含對輸入值的限制。例如，n()必須是整數(shù)，可選參數(shù)必須是實(shí)數(shù)?？蛇x參數(shù)還包括一個(gè)默認(rèn)值，如果用戶在輸入程序名稱時(shí)未*輸入?yún)?shù)，則使用該值。例如，除非用戶*不同的值，否則將為m0分配值0。

注意，在此示例中，n()是樣本大小，alpha()是alpha級別，m0()是假設(shè)零假設(shè)的平均值，ma()是假設(shè)備選假設(shè)的平均值，sd()是標(biāo)準(zhǔn)偏差。

然后drawnorm使用輸入?yún)?shù)從均值為`ma’、標(biāo)準(zhǔn)差為`sd’的正態(tài)分布中生成一個(gè)觀測值為`n’的樣本。接下來，ttest測試樣本均值等于`m0′的原假設(shè)。

可以通過輸入return scalar reject =（r（p）<`alpha'）來返回假設(shè)檢驗(yàn)的結(jié)果。 回想一下ttest將雙邊p值存儲(chǔ)在標(biāo)量r(p)中。 當(dāng)r(p)小于`alpha’*的alpha級別時(shí)，我們的程序使用標(biāo)量拒絕返回值1，否則返回0。

現(xiàn)在，可以輸入simttest和輸入?yún)?shù)來運(yùn)行模擬。







通過輸入return list，可以看到模擬的結(jié)果。







如何使用模擬多次運(yùn)行程序

simttest程序?qū)?zhí)行模擬的一次迭代所需的所有操作。 接下來，需要一種方法來多次運(yùn)行simttest并收集結(jié)果。下面的代碼塊顯示了如何使用simulate來完成這兩項(xiàng)任務(wù)。




參數(shù)reject=r(reject)告訴simulate將r(reject)中返回的結(jié)果保存到名為reject的變量中。 選項(xiàng)reps(100)指示模擬運(yùn)行程序100次。 選項(xiàng)seed(12345)設(shè)置隨機(jī)數(shù)種子，以便我們的結(jié)果可以重現(xiàn)。

冒號(hào)之后是simttest以及我們模擬的輸入?yún)?shù)。 某些模擬需要很長時(shí)間才能運(yùn)行，并且simulate會(huì)在結(jié)果窗口中顯示一個(gè)點(diǎn)，以便知道它仍在運(yùn)行。下面的輸出顯示了我們的模擬結(jié)果。







simulate將結(jié)果保存到變量reject，如果零假設(shè)的測試被拒絕則包含1，否則為0。




可以使用summarize來計(jì)算reject的平均值，該平均值等于100次迭代中拒絕原假設(shè)的次數(shù)所占的比例。根據(jù)輸入?yún)?shù)，該比例是對統(tǒng)計(jì)功效的估算！







在這個(gè)例子中，比例等于0.91，這意味著當(dāng)檢驗(yàn)樣本平均值等于75的替代假設(shè)時(shí)，如果樣本平均值等于70，假設(shè)標(biāo)準(zhǔn)偏差為15且樣本量為100，可以期望91%的功率。

用Power OneMean檢查結(jié)果

可以使用power onemean檢查蒙特卡羅模擬的結(jié)果。 以下示例包含與我們的模擬相同的輸入?yún)?shù)。




power onemean計(jì)算的功率為0.9100，這與我們的模擬估算的功率相同。因?yàn)槟M是隨機(jī)的，所有結(jié)果并不總是**的。改變隨機(jī)數(shù)種子或迭代次數(shù)可以稍微改變估算的功率。但結(jié)果應(yīng)該很接近。

總結(jié)

在本文中，介紹了使用蒙特卡羅模擬計(jì)算統(tǒng)計(jì)功率和樣本大小需求所需的工具。 下一次，將向您展示如何使用Stata的power命令來運(yùn)行模擬，以便可以輕松地為一系列輸入?yún)?shù)創(chuàng)建表格和圖形。








北京天演融智軟件有限公司專注于matlab,stata,spsspro,pscad,gams,nvivo,cyme,gms,mathematica等

• 詞條

  詞條說明

• 2023年P(guān)SCAD軟件用戶大會(huì)視頻回放（一）

  2023年大電網(wǎng)模式下電磁暫態(tài)仿真技術(shù)交流會(huì)暨PSCAD軟件用戶大會(huì)在成都落下帷幕，大會(huì)以“Better Tech Drives The Power World”為主旨，在2天的時(shí)間里為大家?guī)矶鄠€(gè)熱門主題。本系列教程為大會(huì)的回放視頻?？梢悦赓M(fèi)觀看。?首先為大家?guī)淼氖荕anitoba Hydro International Ltd**工程師Shan Jiang和明陽智慧能源集團(tuán)股份公司

• stata軟件免費(fèi)培訓(xùn)

  StataCorp LP公司推出系列免費(fèi)在線培訓(xùn)，較近一期：We have announced a new session of the Ready. Set. Go Stata. webinar. ???? ??? Ready. Set. Go Stata. ??? ??&nb

• NVivo—主題分析比你想象的較受歡迎

  主題分析的實(shí)踐被廣泛應(yīng)用于定性分析，但有時(shí)與其他方法不被承認(rèn)或混淆。在QSR International，我們打破了主題分析作為一種方法的模糊性，并希望這種解釋可以為其注入新活力，因?yàn)樾碌暮驼诔霈F(xiàn)的內(nèi)容形式對已經(jīng)建立的研究工具變得較加不可或缺。什么是主題分析？主題分析不是一種方法，而是一種可以跨不同方法使用的工具（Boyatzis 1998）。 它用于查找內(nèi)容中的常見主題，例如：.文本來源，例如

• CYME Integration Capacity Analysis | 配電系統(tǒng)線路中較大發(fā)電量和電容量分析

  考慮電能質(zhì)量和可靠性，確定配電系統(tǒng)的主機(jī)電容。?在太陽能和風(fēng)力發(fā)電場以及價(jià)格合理的屋頂光伏板等新技術(shù)的支持下，清潔能源的**趨勢改變了配電系統(tǒng)的格局。分布式能源（DER）互連需求的不斷增加催生了對新的創(chuàng)新工程工具的需求。CYME Integration Capacity Analysis模塊可以在不影響配電系統(tǒng)的可靠性和電能質(zhì)量的前提下，對配電系統(tǒng)的發(fā)電或負(fù)載承載能力進(jìn)行評估。