原著者のやる気を削いでしまうかもしれませんが，一般的にプログラムの改良（?）のポイントを提示するという題材とさせてもらいました。

ポイントは

関数化，サブルーチン化する（わかりやすくなる...改良すべき箇所が明瞭になる）
無駄な設定を省く。
余分な変数は使わない（アルゴリズムの改良にもなる）
for ループを避ける（1. にも該当）
その他...細かすぎるので，プログラム参照（Rstudio の reformat code に文句を言わせないようにする）

一枚がビンゴか = function(M, N, Y) {
    # 1枚のビンゴカード X が抽選結果 Y によりビンゴかどうかを判定する
    # ビンゴだったマス目を 0 にし，最後に縦横斜めに 0 が並ぶか判定結果を返す
    X = matrix(sample(M, 25), ncol = 5) # 用紙
    X[3, 3] = 0
    for (i in 1:N) {
        X[X == Y[i]] = 0
    }
    return(0 %in% c(colSums(X), rowSums(X), sum(diag(X)), sum(diag(X[5:1,]))))
    
}

何枚がビンゴか = function(A, M, N, Y) {
    # A 枚のビンゴカードでビンゴになったかどうかを調べ(sub2)，ビンゴだったカード数を返す
    bingo = 0
    for (a in 1:A) {
        bingo = bingo + 一枚がビンゴか(M, N, Y)
    }
    return(bingo)
}

シミュレーション結果の記録 = function(A, B, M, N) {
    # シミュレーションを B 回繰り返し，毎回のビンゴカード枚数を記録する
    # 関数にする（平均値を知りたいなら，平均値を計算して返す方がよい）
	bingos = numeric(B) # シミュレーション結果保存
    Y = sample(M, N) # 抽選番号
	for (b in 1:B) {
		bingos[b] = 何枚がビンゴか(A, M, N, Y) # ビンゴ枚数の記録
	}
	return(bingos)
}

A = 50 # 用紙枚数
B = 5000 # シミュレーション回数
M = 75 # 数字の最大値
N = 25 # 抽選回数

set.seed(123)
system.time({
	print(mean(シミュレーション結果の記録(A, B, M, N)))
})

# 当初プログラム
# [1] 3.2378
#    user  system elapsed 
# 12.008   0.078  12.084 
# 抽選番号は全参加者に対して 1 回のみ（共通） Y が一定になったので，結果は別のものになる
# 3.2084
#    ユーザ   システム       経過  
#     11.062      0.055     11.224 
# 穴あけを X に対して行い，ビンゴ判定も X 自体で行う
# [1] 3.2084
#    ユーザ   システム       経過  
#     10.024      0.032     10.056 
# 要素が 0 を含むかどうかを 0, 1 で返す方法を変える
# [1] 3.2084
#    ユーザ   システム       経過  
#      9.838      0.024      9.861 
# ビンゴのチェックでの二重 for ループをやめる　==> 劇的変化（倍速）
# [1] 3.2084
#    ユーザ   システム       経過  
#      4.742      0.026      4.769 
# 当初の 2.5 倍くらい速くなったかな?
    [1] 3.2084
       ユーザ   システム       経過  
         4.748      0.036      4.785 

system.time({
	print(mean(シミュレーション結果の記録(30, 10000, M, 50)))
})
    [1] 24.4023
       ユーザ   システム       経過  
         7.550      0.034      7.585