「《命》に関わる確率」を疑似乱数を使って解いてみる

package prob

import (
	"math"
	"math/rand"
	"time"
)

func New(p float64) <-chan bool {
	ch := make(chan bool)
	go func() {
		defer close(ch)
		max := 1000000
		limit := percent(p, max)
		rnd := rand.New(rand.NewSource(time.Now().UnixNano()))
		for {
			n := rnd.Intn(max) + 1
			ch <- n < limit
		}
	}()
	return ch
}

func percent(f float64, max int) int {
	if f < 0 {
		return 0
	}
	if f > 1.0 {
		return max
	}
	return int(math.Floor(f*float64(max) + 0.5))
}

ch := prob.New(0.1)

func probability(try int, ch <-chan bool) int {
	count := 0
	for i := 0; i < try; i++ {
		if <-ch {
			count++
		}
	}
	return count
}

func main() {
	ch := prob.New(0.1)
	min := float64(1)
	max := float64(0)
	sum := float64(0)
	sum2 := float64(0)
	try := 10000
	tryf := float64(try)
	ps := make([]float64, 0, try)
	for i := 0; i < try; i++ {
		count := probability(try, ch)
		p := float64(count) / tryf
		ps = append(ps, p)
		if p < min {
			min = p
		}
		if p > max {
			max = p
		}
		sum += p
		sum2 += p * p
	}
	fmt.Printf("minimum value: %7.5f\n", min)
	fmt.Printf("maximum value: %7.5f\n", max)
	ave := sum / tryf
	fmt.Printf("average: %7.5f\n", ave)
	devi := math.Sqrt(sum2/tryf - ave*ave) //standard deviation
	ct := 0
	for _, p := range ps {
		if ave-devi <= p && p <= ave+devi {
			ct++
		}
	}
	fmt.Printf("standard deviation: %7.5f (%4.1f%%)\n", devi, float64(ct)*100.0/tryf)
}

$ go run prob/sample/sample.go
minimum value: 0.08930
maximum value: 0.11160
average: 0.09999
standard deviation: 0.00299 (67.9%)

package note257

import (
	"github.com/spiegel-im-spiegel/mathgirl-problem/prob"
)

type People struct {
	infect <-chan bool
}

type Person bool

func NewPeople() *People {
	return &People{infect: prob.New(0.01)}
}

func (ppl *People) SelectPersion() Person {
	return Person(<-ppl.infect)
}

func (psn Person) Infection() bool {
	return bool(psn)
}

type TestKit struct {
	probability <-chan bool
}

func NewTestKit() *TestKit {
	return &TestKit{probability: prob.New(0.9)}
}

func (tk *TestKit) Inspect(psn Person) bool {
	if psn.Infection() {
		return <-tk.probability
	}
	return !<-tk.probability
}

func probability(ppl *note257.People, tk *note257.TestKit, try int) float64 {
	total := 0
	count := 0
	for i := 0; i < try; i++ {
		psn := ppl.SelectPersion()
		if tk.Inspect(psn) {
			total++
			if psn.Infection() {
				count++
			}
		}
	}
	return float64(count) / float64(total)
}

func main() {
	flag.Parse()
	try, err := strconv.Atoi(flag.Arg(0))
	if err != nil {
		fmt.Fprintln(os.Stderr, err)
		return
	}

	ppl := note257.NewPeople()
	tk := note257.NewTestKit()
	min := float64(1)
	max := float64(0)
	sum := float64(0)
	sum2 := float64(0)
	tryf := float64(try)
	ps := make([]float64, 0, try)
	for i := 0; i < try; i++ {
		p := probability(ppl, tk, try*10)
		ps = append(ps, p)
		if p < min {
			min = p
		}
		if p > max {
			max = p
		}
		sum += p
		sum2 += p * p
	}
	fmt.Printf("minimum value: %7.5f\n", min)
	fmt.Printf("maximum value: %7.5f\n", max)
	ave := sum / tryf
	fmt.Printf("average: %7.5f\n", ave)
	devi := math.Sqrt(sum2/tryf - ave*ave) //standard deviation
	ct := 0
	for _, p := range ps {
		if ave-devi <= p && p <= ave+devi {
			ct++
		}
	}
	fmt.Printf("standard deviation: %7.5f (%4.1f%%)\n", devi, float64(ct)*100.0/tryf)
}

$ go run note257/sample/sample.go 10000

$ go run note257/sample/sample.go 10000
minimum value: 0.08061
maximum value: 0.10279
average: 0.09092
standard deviation: 0.00272 (68.5%)