基于概率模型的蘑菇菌丝规模分析

import math
from random import Random

import numpy as np


def generate_circle_point(r):
    rand = Random()
    angle = rand.random() * (2 * np.pi)
    r_len = r * math.sqrt(rand.random())

    return np.array([r_len * math.sin(angle), r_len * math.cos(angle)])


def generate_circle(m, n, r, point_num):
    points = []
    for i in range(point_num):
        points.append(generate_circle_point(r))

    points = np.array(points)

    points += np.array([m, n])

    return points

import math

import numpy as np


def circle_center1(points):
    return points.mean(axis=0)

def circle_center2(points):
    x_min, y_min = points.min(axis=0)
    x_max, y_max = points.max(axis=0)
    return np.array([(x_min + x_max) / 2, (y_min + y_max) / 2])


def circle_center3(points):
    return np.median(points, axis=0)


def circle_radius1(radius_list):
    r_mean = radius_list.mean()
    return (3 / 2) * r_mean


def circle_radius2(radius_list):
    return max(radius_list)


def circle_radius3(radius_list):
    r_min = min(radius_list)
    r_max = max(radius_list)
    return (r_max + math.sqrt(r_max ** 2 + 4 * r_min ** 2)) / 2

from math_model2 import *

mushrooms = []

with open("mushrooms.txt", "r") as f:
    for i, line in enumerate(f):
        if i == 0:
            continue
        grid = line.split("\t")
        mushrooms.append([float(grid[1]), float(grid[2])])

mushrooms = np.array(mushrooms)

calc_center_funcs = [circle_center1, circle_center2, circle_center3]
calc_radius_funcs = [circle_radius1, circle_radius2, circle_radius3]

with open("ans.txt", "w") as f:
    f.write(f"圆心估计模型\t半径估计模型\t横坐标\t纵坐标\t半径\n")
    for i, center_func in enumerate(calc_center_funcs):
        for j, radius_func in enumerate(calc_radius_funcs):
            m, n = center = center_func(mushrooms)

            r_list = np.linalg.norm(mushrooms - center, axis=1)

            r = radius_func(r_list)

            f.write(f"center{i+1}\tradius{j+1}\t{m:.2f}\t{n:.2f}\t{r:.2f}\n")

import matplotlib.pyplot as plt

from math_model2 import *
from random_generator import generate_circle

num = 20
m = 50
n = 50
r = 20

calc_center_funcs = [circle_center1, circle_center2, circle_center3]
calc_radius_funcs = [circle_radius1, circle_radius2, circle_radius3]

center = [[], [], []]
radius = [[], [], []]

for step in range(200):
    mushrooms = generate_circle(m, n, r, num)
    radius_list = np.linalg.norm(mushrooms - np.array([m, n]), axis=1)

    for i in range(3):
        center_func = calc_center_funcs[i]
        center_ans = center_func(mushrooms)

        center[i].append(center_ans)

    for i in range(3):
        radius_func = calc_radius_funcs[i]
        radius_ans = radius_func(radius_list)

        radius[i].append(radius_ans)

table_val = np.zeros((3, 6))
row_label = ["model 1", "model 2", "model 3"]
col_label = ["x avg", "x error", "x std", "y avg", "y error", "y std"]
for idx, one_center in enumerate(center):
    x, y = zip(*one_center)
    x = np.array(x)
    y = np.array(y)

    table_val[idx][0] = x.mean()
    table_val[idx][1] = x.mean() - m
    table_val[idx][2] = x.std()
    table_val[idx][3] = y.mean()
    table_val[idx][4] = y.mean() - n
    table_val[idx][5] = y.std()

    plt.suptitle(f"Centroid Estimation Model {idx+1}")
    plt.subplot(1, 2, 1)
    plt.hist(x)
    plt.title("Estimation for X Axis")

    plt.subplot(1, 2, 2)
    plt.hist(y)
    plt.title("Estimation for Y Axis")
    plt.show()

table_val = np.around(table_val, decimals=2)
plt.suptitle(f"Numerical Simulation for Centroid Estimation Models")
plt.table(
    cellText=table_val,
    rowLabels=row_label,
    colLabels=col_label,
    loc="center",
    cellLoc="center",
    rowLoc="center",
)
plt.axis("off")
plt.show()

table_val = np.zeros((3, 3))
row_label = ["model 1", "model 2", "model 3"]
col_label = ["r avg", "r error", "r std"]
for idx, one_radius in enumerate(radius):
    one_radius = np.array(one_radius)

    table_val[idx][0] = one_radius.mean()
    table_val[idx][1] = one_radius.mean() - r
    table_val[idx][2] = one_radius.std()
    
    plt.hist(one_radius)
    plt.suptitle(f"Radius Estimation Model {idx+1}")
    plt.show()

table_val = np.around(table_val, decimals=2)
plt.suptitle("Numerical Simulation for Radius Estimation Models")
plt.table(
    cellText=table_val,
    rowLabels=row_label,
    colLabels=col_label,
    loc="center",
    cellLoc="center",
    rowLoc="center",
)
plt.axis("off")
plt.show()