Forum: Ders Arası RSS
Monte Carlo tümlevi
erdem (Moderatör) #1
Üye Tem 2009 tarihinden beri · 978 mesaj · Konum: Eskişehir
Grup üyelikleri: Genel Moderatörler, Üyeler
Profili göster · Bu konuya bağlantı
Konu adı: Monte Carlo tümlevi
Dün konuştuğumuz Monte Carlo tümlev alma yöntemi hakkında birisi bir Python programı yazmış. Sonra Riemann yönteminin daha hızlı ve daha doğru sonuçlar verdiğinden bahsetmiş.

Ben de hızlıca bir D kodu yazayım dedim.

import std.stdio, std.random;
 
void main()
{
    int icerdekiNoktalar = 0, toplamNokta = 0;
 
    Random cekirdek;
    for (int i = 0; i < 100; ++i)
    {
        auto x = uniform(0.0L, 1.0L, cekirdek);
        auto y = uniform(0.0L, 1.0L, cekirdek);
 
        ++toplamNokta;
        writeln("toplam nokta", toplamNokta);
 
        writeln("x,y ", x," ", y);
 
 
        if (((x * x) + (y * y)) <= 1)
        {
            icerdekiNoktalar++;
            writeln(icerdekiNoktalar);
        }
    }
 
    writeln("icerdeki ", icerdekiNoktalar);
    writeln("toplam ", toplamNokta);
 
    float pisayisi = cast (float) (icerdekiNoktalar/toplamNokta);
 
 
    writeln("Pi sayisi = ", pisayisi);
 
}

Ama nedense pi sayısının çeyreğini bulurken sonuç 0 çıkıyor.

Sanırım biraz da D çalışmanın zamanı gelmiş ;)
Avatar
Salih Dinçer #2
Üye Ock 2012 tarihinden beri · 1912 mesaj · Konum: İstanbul
Grup üyelikleri: Üyeler
Profili göster · Bu konuya bağlantı
Sanırım şu şekilde yapmalısın:
    float pisayisi = cast(float)icerdekiNoktalar/cast(float)toplamNokta;
Ama belki de hiç dönüştürmeye gerek yoktur. Bütün değişkenler double olamaz mı?
Bilgi paylaştıkça bir bakmışız; kar topu olmuş ve çığ gibi üzerimize geliyor...:)
Avatar
Salih Dinçer #3
Üye Ock 2012 tarihinden beri · 1912 mesaj · Konum: İstanbul
Grup üyelikleri: Üyeler
Profili göster · Bu konuya bağlantı
Oluyormuş:
/*
 * mcPi.d: Monte Carlo Pi
 * Source: http://ddili.org/forum/thread/1292
 * Programmer: Erdem Öncel
 */
import std.stdio, std.random;
 
void main()
{
    double icerdekiNoktalar = 0.0, toplamNokta = 0.0;
    double resolution = 2000.0;
    
    Random cekirdek;
    do {
        auto x = uniform(0.0L, 1.0L, cekirdek);
        auto y = uniform(0.0L, 1.0L, cekirdek);
 
        toplamNokta++;
        debug(1)
        {
            writeln("toplam nokta", toplamNokta);
            writeln("x, y ", x," ", y);
        }
 
        if (((x * x) + (y * y)) <= 1)
        {
            icerdekiNoktalar++;
            debug(2) writeln(icerdekiNoktalar);
        }
    } while(icerdekiNoktalar < resolution);
 
    writefln("icerdeki: %.0f", icerdekiNoktalar);
    writefln("toplam: %.0f", toplamNokta);
 
    auto pisayisi = (icerdekiNoktalar/toplamNokta) * 4.0;
 
    writeln("Pi sayisi = ", pisayisi);
}
Bilgi paylaştıkça bir bakmışız; kar topu olmuş ve çığ gibi üzerimize geliyor...:)
erdem (Moderatör) #4
Üye Tem 2009 tarihinden beri · 978 mesaj · Konum: Eskişehir
Grup üyelikleri: Genel Moderatörler, Üyeler
Profili göster · Bu konuya bağlantı
Teşekkürler double kullanınca çalıştı.

Ama üçüncü basamağa kadar doğru tahmin edebilmesi için bir milyondan fazla tahmin yapması gerekiyor.

Hatta aslında sorduğum soru buydu. Birisi videoyla cevap vermiş.

http://www.youtube.com/watch?v=djJzd-MtVgM

Orada ilk anlattığı benim çözümüm ;-) Ama ikinci çözüm de güzel.
Bu mesaj erdem tarafından değiştirildi; zaman: 2013-11-08, 12:28.
Avatar
Salih Dinçer #5
Üye Ock 2012 tarihinden beri · 1912 mesaj · Konum: İstanbul
Grup üyelikleri: Üyeler
Profili göster · Bu konuya bağlantı
Bende bir kaç milyonu denedim ve hislerim Pi'den uzaklaştığını söylüyor...

Belki de beşik gibi bir salınımı olabilir. Çünkü bir ara neredeyse 3.14'ün altına indiğini ve sonra yükseldiğini gördüm. Belki belli aralıklarla numune alınıp (örneğin 100 bin'de 1) değerler bir çizelgeye aktarılabilir. Çıkan eğri çok ilginç olabilir...
Bilgi paylaştıkça bir bakmışız; kar topu olmuş ve çığ gibi üzerimize geliyor...:)
acehreli (Moderatör) #6
Kullanıcı başlığı: Ali Çehreli
Üye Haz 2009 tarihinden beri · 4527 mesaj
Grup üyelikleri: Genel Moderatörler, Üyeler
Profili göster · Bu konuya bağlantı
Bir tane de yalnızca Phobos ve aralıklarla:
import std.stdio;
import std.range;
import std.algorithm;
import std.typecons;
import std.random;
import std.conv;
 
void main()
{
    enum size_t toplam = 1_000_000;
 
    const size_t içerdekiler =
        iota(toplam)
        .map!(_ => tuple(uniform(0.0L, 1.0L), uniform(0.0L, 1.0L)))
        .count!(a => ((a[0] * a[0]) + (a[1] * a[1])) < 1);
 
    const pi = 4.0 * içerdekiler / toplam;
 
    /* Başka yollarla da olur. Bir de std.conv.to ile:
     *
     *   const pi = içerdekiler.to!double / toplam * 4;
     */
 
    writeln(pi);
}
Ali
erdem (Moderatör) #7
Üye Tem 2009 tarihinden beri · 978 mesaj · Konum: Eskişehir
Grup üyelikleri: Genel Moderatörler, Üyeler
Profili göster · Bu konuya bağlantı
Çok güzel program. Teşekkürler! :)

Ali beyin yazdığı programı 10 milyon sayıyla denedim üçüncü ondalık basamağa kadar buldu.

Python programını yazan kişininkini de denedim. O da üçüncü ondalık basamağa kadar buldu.

Pi sayısı 3.14159 şeklinde devam ediyor. Gerçekten de Riemann yöntemi daha hızlı çalışıyor ama daha hassas bir sonuç bulamıyor.

Merak edenler için Python programını da yazayım.

import math
import random
 
# from time import perf_counter
# burayı alttaki satırla değiştirdim
 
from timeit import default_timer as perf_counter
 
def monte_carlo(start, stop, num_points, f):
    hits = 0
    upp_bound = 0
    for i in range(num_points):
        x = random.uniform(start, stop)
        y = random.uniform(0, 4) # this is cheating since I already know it's 4
        if y <= f(x):
            hits += 1
        if y > upp_bound:
            upp_bound = y
    ans = (hits / num_points) * ((stop - start) * upp_bound)
    return ans
 
def riemann(start, stop, step, f):
    return sum(f(start+step*m)*step for m in range(int((stop - start)/step)))
 
f = lambda x: math.sqrt(16 - 16*x**2)
 
start = perf_counter()
a = monte_carlo(0, 1, 10000000, f)
m_time = perf_counter() - start
start = perf_counter()
b = riemann(0, 1, 1/15000, f)
r_time = perf_counter() - start
print('Monte Carlo: %s\nTime elapsed: %s' % (a, m_time))
print('Riemann sum: %s\nTime elapsed: %s' % (b, r_time))

Burada işleve bir fonksiyon göndererek direkt olarak bahsettiğim sorunun integralini hesaplatmış anladığım kadarıyla.
Avatar
Salih Dinçer #8
Üye Ock 2012 tarihinden beri · 1912 mesaj · Konum: İstanbul
Grup üyelikleri: Üyeler
Profili göster · Bu konuya bağlantı
Bende istatistiksel açıdan inceleyim dedim...:)

Malum hissettiğim bir düzensizlik vardı. Bunun matematiksel karşılığı nedir bilmiyorum ama bazen +/- dalgalanma (peak) yapıyor. Tabi sonuçları, örnekleme çözünürlüğü on binde bir olacak şekilde çizelgeye yansıttığım için derin/yüksek tepeler pek belli olmuyor. Kullandığım kod ise aşağıda ve bunu döngü içindeki koşula (denklem kapsamına) koymalısınız...
            enum çözünürlük = 10_000.0;
            if(içerdekiNoktalar % çözünürlük == 0) {
              auto örnekleme = (3.141592 * çözünürlük) -
                               (içerdekiNoktalar/toplamNokta * 4.0 * çözünürlük);
              //if(örnekleme > - çözünürlük && örnekleme < çözünürlük)
              // arada dalgalanma (peak) var!/neden düzensiz olabilir ki?
              writefln("%.0f", örnekleme);
             
            }
[Resim: http://img27.imageshack.us/img27/1205/lcv3.png]

Bu çizelgeyi, elimden geldiğince yumuşatmaya ve gereksiz ayrıntılardan uzak tutmaya çalıştım. Gri kutu ile işaretlediğim başlangıçtaki kararsızlık normal çünkü henüz iç/dış noktaların oranı yeterince büyük sayı değil. Ama sonradan 1 milyon döngüye rağmen, bazen Pi sayısına yaklaşıyor, bazense uzaklaşıyor!

Dip Not: Tam olarak denk geldiği noktayı kestirmek güç olduğu için başka algoritmalar ile desteklemeden Pi sayısını bulmak için kullanılabilir gözükmüyor. Ama rasgeleliğin sonuçlarını değerlendirmek açısından güzel...:)

Sevgiler, saygılar...
Bilgi paylaştıkça bir bakmışız; kar topu olmuş ve çığ gibi üzerimize geliyor...:)
Doğrulama Kodu: VeriCode Lütfen resimde gördüğünüz doğrulama kodunu girin:
İfadeler: :-) ;-) :-D :-p :blush: :cool: :rolleyes: :huh: :-/ <_< :-( :'( :#: :scared: 8-( :nuts: :-O
Özel Karakterler:
Forum: Ders Arası RSS
Bağlı değilsiniz. · Şifremi unuttum · ÜYELİK
This board is powered by the Unclassified NewsBoard software, 20100516-dev, © 2003-10 by Yves Goergen
Şu an: 2017-11-18, 22:50:29 (UTC -08:00)